Loading required packages
#importing libraries
library("tidyverse")
library(here)
library(ggrepel)
library(lme4)
library(lmerTest)
library(sjstats)
library(emmeans)
library(MuMIn)
library(sjPlot)
library(sjmisc)
library("car")
Reading data into data frame
#Ethnicity data set
ethnicity_data <- read_delim(file = here("data", "qry_Ethnic_origin_2.txt"), delim = ",")
Rows: 44191 Columns: 3
── Column specification ───────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (2): GENDER_DESC, ETHNIC_ORIGIN_DESC
dbl (1): ID2
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(ethnicity_data)
#Mark data set
mark_data <- read_delim(file = here("data", "qry_Mark_Data_2.txt"), delim = ",")
Rows: 1117342 Columns: 22
── Column specification ───────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (13): PROGRAMME_SCHOOL, PROGRAMME_OWNER, PROGRAMME_LEVEL_DESC, PROGRAMME_PART, PART_END...
dbl (9): ID2, STUDENT_START_YEAR, PART_START_ACADEMIC_YEAR, INSTANCE_ACADEMIC_YEAR, MODULE...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(mark_data)
#checking for missing data
ethnicity_data %>% sapply(function(x) sum(is.na(x)) )
ID2 GENDER_DESC ETHNIC_ORIGIN_DESC
0 1 15
#removing missing data
clean_ethnicity_data <- ethnicity_data %>%
drop_na()
#checking if na's removed.
clean_ethnicity_data %>% sapply(function(x) sum(is.na(x)))
ID2 GENDER_DESC ETHNIC_ORIGIN_DESC
0 0 0
#checking for duplicates & removing
length(unique(clean_ethnicity_data$ID2)) #number of students in ethnic data
[1] 44176
length(unique(clean_ethnicity_data$ID2)) == nrow(clean_ethnicity_data) #returns TRUE as no duplicates
[1] TRUE
#renaming columns
clean_ethnicity_data %>%
rename(
Gender = GENDER_DESC,
Ethnicity = ETHNIC_ORIGIN_DESC
) -> clean_ethnicity_data
#count of ethnic groups / adding column for % of each ethnicity
#check levels of gender
clean_ethnicity_data %>%
group_by(Gender) %>%
tally () %>%
mutate(freq = (n / sum(n))*100) %>%
mutate(percentage = scales::label_percent()(n / sum(n))) %>%
rename("Number of students" = n)
#Categorising into man vs woman other levels to little cases
#remove non binary, other & prefer not to say
clean_ethnicity_data<-clean_ethnicity_data[!(clean_ethnicity_data$Gender=="Non-Binary" | clean_ethnicity_data$Gender=="Prefer not to say" | clean_ethnicity_data$Gender=="Other") ,]
#check levels of ethnicites
clean_ethnicity_data %>%
group_by(Ethnicity) %>%
tally () %>%
mutate(freq = (n / sum(n))*100) %>%
mutate(percentage = scales::label_percent()(n / sum(n)))
I will categorise ethnicity differently depending on research question. Firstly, I will begin with white vs other.
#duplicating data set at this point as will categorise differently for second research question
clean_ethnicity_data -> clean_ethnicity_data2
#categorise ethnicities into white vs non white
#rename all white columns- white.
#rename all other columns other.
clean_ethnicity_data[clean_ethnicity_data == "White - British"] <- "White"
clean_ethnicity_data[clean_ethnicity_data == "White - Irish"] <- "White"
clean_ethnicity_data[clean_ethnicity_data == "White - Scottish"] <- "White"
clean_ethnicity_data[clean_ethnicity_data == "Other White Background"] <- "White"
clean_ethnicity_data[clean_ethnicity_data == "Arab"] <- "BAME"
clean_ethnicity_data[clean_ethnicity_data == "Asian Or Asian British - Bangladeshi"] <- "BAME"
clean_ethnicity_data[clean_ethnicity_data == "Asian Or Asian British - Indian"] <- "BAME"
clean_ethnicity_data[clean_ethnicity_data == "Asian Or Asian British - Pakistani"] <- "BAME"
clean_ethnicity_data[clean_ethnicity_data == "Black Or Black British - African"] <- "BAME"
clean_ethnicity_data[clean_ethnicity_data == "Black Or Black British - Caribbean"] <- "BAME"
clean_ethnicity_data[clean_ethnicity_data == "Chinese"] <- "BAME"
clean_ethnicity_data[clean_ethnicity_data == "Gypsy or Traveller"] <- "BAME"
clean_ethnicity_data[clean_ethnicity_data == "Mixed - White And Asian"] <- "BAME"
clean_ethnicity_data[clean_ethnicity_data == "Mixed - White And Black African"] <- "BAME"
clean_ethnicity_data[clean_ethnicity_data == "Mixed - White And Black Caribbean"] <- "BAME"
clean_ethnicity_data[clean_ethnicity_data == "Other Asian Background"] <- "BAME"
clean_ethnicity_data[clean_ethnicity_data == "Other Black Background"] <- "BAME"
clean_ethnicity_data[clean_ethnicity_data == "Other Ethnic Background"] <- "BAME"
clean_ethnicity_data[clean_ethnicity_data == "Other Mixed Background"] <- "BAME"
#remove not known,no information & info refused- as could be any ethnicity.
clean_ethnicity_data<-clean_ethnicity_data[!(clean_ethnicity_data$Ethnicity=="Information Refused" | clean_ethnicity_data$Ethnicity=="Information Refused - OBSOLETE VALUE" | clean_ethnicity_data$Ethnicity=="No Information" | clean_ethnicity_data$Ethnicity=="Not Known") ,]
#changing numerical variables to factors.
clean_ethnicity_data %>%
mutate_at(vars(Gender, Ethnicity), list(factor)) -> clean_ethnicity_data
summary(clean_ethnicity_data)
ID2 Gender Ethnicity
Min. : 1 Man :24727 BAME :17389
1st Qu.:11045 Woman:18025 White:25363
Median :22124
Mean :22101
3rd Qu.:33137
Max. :44191
#Ethnicity
#Visualising number of students in each ethnicity category
ethnicity_data %>%
group_by(ETHNIC_ORIGIN_DESC) %>%
summarise(counts = n()) %>%
mutate(freq = (counts / sum(counts))*100) %>%
mutate(percentage = scales::label_percent()(counts / sum(counts))) %>%
ggplot(., aes(x = ETHNIC_ORIGIN_DESC, y = counts)) +
geom_bar(fill = "#0073C2FF", stat = "identity") +
#geom_text(aes(label = percentage), angle = 60) +
theme(axis.text.x = element_text(angle = 60,
hjust = 1)) +
ggtitle("Number of students in each ethnic category") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab('Ethnicity') +
ylab('Number of students')
#Ethnicity
#Percentage of students belonging to each ethnicity after categorisation
clean_ethnicity_data %>%
group_by(Ethnicity) %>%
tally () %>%
mutate(freq = scales::label_percent()(n / sum(n))) %>%
ggplot( ., aes(x="", y=freq, fill=Ethnicity)) +
geom_bar(stat="identity", width=1, color="white") +
coord_polar("y", start=0) +
theme_void() +
geom_text(aes(label = freq), color = "white", position = position_stack(vjust = 0.5), size = 4) +
scale_fill_brewer(palette = "Paired") +
ggtitle("Percentage of students in each ethnic category") +
theme(plot.title = element_text(hjust = 0.5))
#GENDER
#Percentage of students belonging to each gender category
ethnicity_data %>%
group_by(GENDER_DESC) %>%
summarise(counts = n()) %>%
mutate(freq = (counts / sum(counts))*100) %>%
mutate(percentage = scales::label_percent()(counts / sum(counts))) %>%
ggplot(., aes(x = GENDER_DESC, y = counts, fill=GENDER_DESC)) +
geom_bar(stat = "identity") +
geom_text(aes(label = percentage), vjust = -0.5) +
scale_fill_brewer(palette = "Paired") +
ggtitle("Percentage of students in each gender category") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab('Gender') +
ylab('Number of students') +
guides(fill = guide_legend(title = "Gender Description "))
#GENDER
#Percentage of students belonging to each gender after categorisation
clean_ethnicity_data %>%
group_by(Gender) %>%
tally () %>%
mutate(freq = scales::label_percent()(n / sum(n))) %>%
ggplot( ., aes(x="", y=freq, fill=Gender)) +
geom_bar(stat="identity", width=1, color="white") +
coord_polar("y", start=0) +
theme_void() +
geom_text(aes(label = freq), color = "white", position = position_stack(vjust = 0.5), size = 5) +
scale_fill_brewer(palette = "Set2") +
ggtitle("Percentage of students in each gender category") +
theme(plot.title = element_text(hjust = 0.5))
#summary of mark data
summary(mark_data)
ID2 STUDENT_START_YEAR PROGRAMME_SCHOOL PROGRAMME_OWNER PROGRAMME_LEVEL_DESC
Min. : 1 Min. : 7.0 Length:1117342 Length:1117342 Length:1117342
1st Qu.:11013 1st Qu.:17.0 Class :character Class :character Class :character
Median :22121 Median :18.0 Mode :character Mode :character Mode :character
Mean :22094 Mean :17.8
3rd Qu.:33118 3rd Qu.:19.0
Max. :44191 Max. :22.0
PROGRAMME_PART PART_START_ACADEMIC_YEAR PART_END_DATE MODULE_OWNER
Length:1117342 Min. : 8.00 Length:1117342 Length:1117342
Class :character 1st Qu.:18.00 Class :character Class :character
Mode :character Median :19.00 Mode :character Mode :character
Mean :18.66
3rd Qu.:20.00
Max. :21.00
MODULE_SCHOOL MODULE_CODE MODULE_LONG_TITLE INSTANCE_ACADEMIC_YEAR
Length:1117342 Length:1117342 Length:1117342 Min. :17.0
Class :character Class :character Class :character 1st Qu.:18.0
Mode :character Mode :character Mode :character Median :19.0
Mean :18.7
3rd Qu.:20.0
Max. :21.0
MODULE_INSTANCE_NO MODULE_MARK INCLUDE_MARK_DESC MODULE_CREDITS
Min. : 1.000 Min. : 0.00 Length:1117342 Min. : 10.0
1st Qu.: 1.000 1st Qu.: 58.00 Class :character 1st Qu.: 10.0
Median : 1.000 Median : 65.00 Mode :character Median : 15.0
Mean : 1.217 Mean : 63.94 Mean : 16.4
3rd Qu.: 1.000 3rd Qu.: 72.00 3rd Qu.: 20.0
Max. :59.000 Max. :100.00 Max. :120.0
NA's :8
ASSESSMENT_COMPONENT_CODE ASSESSMENT_COMPONENT_TITLE ASSESSMENT_COMPONENT_TYPE_DESC
Length:1117342 Length:1117342 Length:1117342
Class :character Class :character Class :character
Mode :character Mode :character Mode :character
ASSESSMENT_WEIGHT MARK
Min. : 0.00 Min. : -3.00
1st Qu.: 20.00 1st Qu.: 58.00
Median : 40.00 Median : 66.00
Mean : 48.34 Mean : 66.18
3rd Qu.: 75.00 3rd Qu.: 75.00
Max. :100.00 Max. :100.00
#checking for missing data?
mark_data %>% sapply(function(x) sum(is.na(x)) )
ID2 STUDENT_START_YEAR PROGRAMME_SCHOOL
0 0 20509
PROGRAMME_OWNER PROGRAMME_LEVEL_DESC PROGRAMME_PART
0 0 0
PART_START_ACADEMIC_YEAR PART_END_DATE MODULE_OWNER
0 130170 0
MODULE_SCHOOL MODULE_CODE MODULE_LONG_TITLE
20509 0 0
INSTANCE_ACADEMIC_YEAR MODULE_INSTANCE_NO MODULE_MARK
0 0 8
INCLUDE_MARK_DESC MODULE_CREDITS ASSESSMENT_COMPONENT_CODE
0 0 0
ASSESSMENT_COMPONENT_TITLE ASSESSMENT_COMPONENT_TYPE_DESC ASSESSMENT_WEIGHT
0 0 0
MARK
0
#remove missing data??? - decided to not remove missing data as only missing data in programme school, module school and part_end_date and module mark- i am not directly assessing these columns- so doesn't matter too much. If looking at effect of school on mark may need to consider removing missing values.
#checking for duplicates
length(unique(mark_data$ID2)) #number of students in marks data- same as ethnicity data set
[1] 44191
length(unique(mark_data$ID2)) == nrow(mark_data) #returns FALSE as there are duplicates
[1] FALSE
#removing duplicate rows- 1,114,911 rows once duplicates removed.
distinct(mark_data) -> clean_mark_data
#changing variables to factors.
clean_mark_data %>%
mutate_at(vars(PROGRAMME_SCHOOL, PROGRAMME_OWNER, PROGRAMME_LEVEL_DESC, PROGRAMME_PART, PART_END_DATE, MODULE_OWNER, MODULE_SCHOOL, MODULE_CODE, MODULE_LONG_TITLE, INCLUDE_MARK_DESC, ASSESSMENT_COMPONENT_CODE, ASSESSMENT_COMPONENT_TITLE, ASSESSMENT_COMPONENT_TYPE_DESC), list(factor)) -> clean_mark_data
#how many years of data?
clean_mark_data %>%
group_by(STUDENT_START_YEAR) %>%
count ()
#data spans 2007-2022
#Only interested in the last 5 years, as before 2016, a lot less data/outdated.
#mark data only has a single students results for 2022, therefore, only one gender, one ethnicity.
#This is not good enough for any analysis on students performance, will cause serious over generalisation.
#remove data before 2017 & after 2021
clean_mark_data %>%
filter(STUDENT_START_YEAR >= 17 & STUDENT_START_YEAR <= 21) ->clean_mark_data
#Visualisation of students in each year group.
clean_mark_data %>%
group_by(STUDENT_START_YEAR) %>%
count() %>%
ggplot(., aes(x = STUDENT_START_YEAR, y= n)) +
geom_bar(stat = "identity", fill = "darkgreen", size = 2.5) +
geom_text(aes(label = n), vjust = -0.5) +
ggtitle("Number of students in each year (2017-2021)") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab('Student Start Year') +
ylab('Number of students')
#now only 886,263 rows
#How many assessment types in data set?
clean_mark_data %>%
group_by(ASSESSMENT_COMPONENT_TYPE_DESC) %>%
count()
#Coursework, Exam, In-Class Test, In-Dept Exam, Laboratory Test
#types of assessment types & relative percentages of students who did the AT
clean_mark_data %>%
group_by(ASSESSMENT_COMPONENT_TYPE_DESC) %>%
summarise(counts = n()) %>%
mutate(freq = (counts / sum(counts))) %>%
mutate(percentage = scales::label_percent()(counts / sum(counts))) %>%
ggplot( ., aes(x="", y=freq, fill=ASSESSMENT_COMPONENT_TYPE_DESC)) +
geom_bar(stat="identity", width=1, color="white") +
coord_polar("y", start=180) +
theme_void() +
geom_text(aes(label = percentage), position = position_stack(vjust = 0.5), size = 3) +
scale_fill_brewer(palette = "Spectral") +
guides(fill = guide_legend(title = "Assessment type")) +
ggtitle("Percentage of students in each assessment type") +
theme(plot.title = element_text(hjust = 0.5))
#fix labels?? or change to bar chart
#change to dark2?
#categorise into exam and cw only.
#So- exam = Exam, In-Dept Exam, cw = CW, In-Class Test and Lab Test
clean_mark_data[clean_mark_data == "In-Dept Exam"] <- "Exam"
Warning: stack imbalance in '==', 41 then 43
Warning: stack imbalance in '==', 25 then 27
Warning: stack imbalance in '[<-', 10 then 12
Warning: stack imbalance in '<-', 2 then 4
clean_mark_data[clean_mark_data == "In-Class Test"] <- "Coursework"
clean_mark_data[clean_mark_data == "Laboratory Test"] <- "Coursework"
#dropping levels
#in dept exam, in class test and lab test
droplevels(clean_mark_data$ASSESSMENT_COMPONENT_TYPE_DESC) -> clean_mark_data$ASSESSMENT_COMPONENT_TYPE_DESC
summary(clean_mark_data$ASSESSMENT_COMPONENT_TYPE_DESC)
Coursework Exam
618851 267412
#visualising percentage of students under each assessment type after categorisation
clean_mark_data %>%
group_by(ASSESSMENT_COMPONENT_TYPE_DESC) %>%
summarise(counts = n()) %>%
mutate(freq = (counts / sum(counts))) %>%
mutate(percentage = scales::label_percent()(counts / sum(counts))) %>%
ggplot( ., aes(x="", y=freq, fill=ASSESSMENT_COMPONENT_TYPE_DESC)) +
geom_bar(stat="identity", width=1, color="white") +
coord_polar("y", start=0) +
theme_void() +
geom_text(aes(label = percentage), position = position_stack(vjust = 0.5), size = 3) +
scale_fill_brewer(palette = "Spectral") +
guides(fill = guide_legend(title = "Assessment type")) +
ggtitle("Percentage of students in each assessment type") +
theme(plot.title = element_text(hjust = 0.5))
#cleaning programme part
summary(clean_mark_data$PROGRAMME_PART)
A B C D F T
343429 229228 97993 6310 18717 174344
ggplot(clean_mark_data, aes(x = PROGRAMME_PART)) +
geom_bar(fill = "blue") +
theme(axis.text.x = element_text(angle = 60,
hjust = 1)) +
ggtitle("Number of students in each part") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab('Programme Part') +
ylab('Number of students')
As can be seen from the above table- there are many parts- but in this study i am only interested in UG and PG students, so will remove all cases where not A, B, C, D, F, and T.
#as can be seen from the above table- there are many parts- but in this study i am only interested in UG and PG students, so will remove all cases where not A, B, C, D, F, and T.
#There are 16 parts- want to keep 6 and remove 10
subset(clean_mark_data, PROGRAMME_PART != "G" & PROGRAMME_PART != "I" & PROGRAMME_PART != "R" & PROGRAMME_PART != "R0" & PROGRAMME_PART != "R1" & PROGRAMME_PART != "R2" & PROGRAMME_PART != "RT" & PROGRAMME_PART != "TO" & PROGRAMME_PART != "TR" & PROGRAMME_PART != "TX") -> clean_mark_data
#dropping levels with 0
summary(clean_mark_data$PROGRAMME_PART)
A B C D F G I R R0 R1 R2 RT T
343430 229228 97993 6310 18717 0 0 0 0 0 0 0 174344
TO TR TX
0 0 0
droplevels(clean_mark_data$PROGRAMME_PART) -> clean_mark_data$PROGRAMME_PART
summary(clean_mark_data$PROGRAMME_PART)
A B C D F T
343430 229228 97993 6310 18717 174344
ggplot(clean_mark_data, aes(x = PROGRAMME_PART)) +
geom_bar(fill = "blue") +
theme(axis.text.x = element_text(angle = 60,
hjust = 1)) +
ggtitle("Number of students in each part") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab('Programme Part') +
ylab('Number of students')
#distribution of marks in cleaned mark dataset
hist(clean_mark_data$MARK,
main="Histogram of mark (2017-2021)",
xlab="Mark",
xlim=c(0,100),
col="darkmagenta")
#seems close to normal distribution
#from above looks like there are marks below 0
clean_mark_data %>%
group_by(MARK) %>%
tally ()
#REMOVE MARKS BELOW, ONE ENTRY OF -3.
filter(clean_mark_data, MARK > 0) -> clean_mark_data
#now 870021 rows
#JOIN ETHNIC DATA AND MARK DATA
inner_join(clean_ethnicity_data, clean_mark_data) -> full_data
Joining, by = "ID2"
print(distinct(full_data))
#840368 rows
This is a multi-leveled question. I will began by looking at the effect of ethnicity and assessment type separately on student performance.
Variance in marks for white vs other
boxplot(MARK~Ethnicity, data = full_data)
looks like a lot of data points are outliers- but data is over 5 years and we are mostly interested in top grades so will filter these out later.
Analysis of average Marks of White vs BAME/Other students (2017-2021)
Below is bar chart showing of effect of ethnicity on mean mark between 2017-2021
full_data %>%
select(Ethnicity, MARK) %>%
group_by(Ethnicity) %>%
tally(mean(MARK)) %>%
rename(Average_mark = n) %>%
ggplot(aes(x=Ethnicity, y = Average_mark, fill =Ethnicity)) +
geom_bar(stat = 'identity', position = 'dodge') +
scale_fill_brewer(palette = "Paired") +
geom_text(aes(label = round(Average_mark, 2)), position=position_dodge(width=0.9), vjust=-0.25) +
ggtitle("Average Test Scores by Ethnicity between 2017-2021") +
theme(plot.title = element_text(hjust = 0.5)) +
ylab('Average Mark')
Next, I wanted to see the affect of ethnicity on average mark each year between 2017-2021. This table below shows the gap in average grade in white vs Bame students is reducing over the years. This reduction is very small changing by just 1.4 over 5 years.
#White vs non white grades table like in screenshot in report
full_data %>%
select(INSTANCE_ACADEMIC_YEAR, Ethnicity,MARK) %>%
group_by(INSTANCE_ACADEMIC_YEAR, Ethnicity) %>%
tally(mean(MARK)) %>%
rename(Average_mark = n) %>%
spread(key = Ethnicity, value = Average_mark) %>%
mutate(Gap = White - BAME)
Below is a line graph to visualise the above table. This is showing gap is reducing, but both white and bame students average grade seems to be reducing over the years.
full_data %>%
select(INSTANCE_ACADEMIC_YEAR, Ethnicity,MARK) %>%
group_by(INSTANCE_ACADEMIC_YEAR, Ethnicity) %>%
tally(mean(MARK)) %>%
rename(Average_mark = n) %>%
#spread(key = Ethnicity, value = Average_mark) %>%
# mutate(Gap = White - BAME) %>%
ggplot(aes(x = INSTANCE_ACADEMIC_YEAR, y = Average_mark, col = Ethnicity)) +
geom_line() +
geom_point() +
ggtitle("Average grade for white vs BAME students (2017-2021)") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab('Academic Year') +
ylab('Average Mark') +
labs(col='Ethnicity')
The below bar graph is showing percentage of white vs BAME students achieving top degree (2017-2021) NOTE- for methodology/results I grouped by ID2 first in order to calculate average grade for each student..then i filtered for students who got an average grade above 60 Then grouped by ethnicity This is important as the data is not important, if not grouped by ID2 results would show any student who got above 60 for at least 1 assessment
full_data %>%
group_by(ID2) %>%
mutate(Average_mark = mean(MARK)) %>%
filter(Average_mark >=60) %>%
select(ID2, Ethnicity, Average_mark) %>%
distinct() %>%
group_by(Ethnicity) %>%
summarise(counts = n()) %>%
mutate(percentage = scales::label_percent()(counts / sum(counts))) %>%
ggplot(aes(x=Ethnicity, y = counts, fill =Ethnicity)) +
geom_bar(stat = 'identity', position = 'dodge') +
scale_fill_brewer(palette = "Pastel1") +
geom_text(aes(label = percentage), position = position_stack(vjust = 1.1), size = 3) +
ggtitle("Percentage of white vs BAME students achieving top degree (2017-2021)") +
theme(plot.title = element_text(hjust = 0.5)) +
ylab('Number of students')
Below graph shows difference between number of White vs BAME students achieving top degree (BAME awarding gap) is consistent over years 2017-2021
#Visualising percentage of white vs BAME students achieving top degree each year (2017-2021)
full_data %>%
group_by(ID2) %>%
mutate(Average_mark = mean(MARK)) %>%
filter(Average_mark >=60) %>%
select(ID2, Ethnicity, Average_mark, INSTANCE_ACADEMIC_YEAR) %>%
distinct() %>%
group_by(INSTANCE_ACADEMIC_YEAR, Ethnicity) %>%
summarise(counts = n()) %>%
mutate(percentage = scales::label_percent()(counts / sum(counts))) %>%
ggplot(aes(x=INSTANCE_ACADEMIC_YEAR, y = counts, fill =Ethnicity)) +
geom_bar(stat = 'identity', position = 'dodge') +
scale_fill_brewer(palette = "Pastel1") +
geom_text(aes(label = percentage), position=position_dodge(width=0.9), vjust=-0.25) +
ggtitle("Percentage of white vs BAME students achieving top degree each year (2017-2021)") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab('Academic Year') +
ylab('Number of students')
`summarise()` has grouped output by 'INSTANCE_ACADEMIC_YEAR'. You can override using the `.groups` argument.
#this shows awarding gap is reducing but only by 4% between 2017-2020, 2021 less data so not good for comparison
Problem with above graphs is that there are a lot more white students in the university, and therefore it is a given that more white students are achieving top degrees. So instead, this graph clearly shows distribution of degree outcomes instead.The following graph makes the results clearer for interpretation/comparisons. For example here, 37% of white students get a first in comparison to only 21.9% of BAME students.
full_data %>%
group_by(ID2) %>%
mutate(Average_mark = mean(MARK)) %>%
mutate(Grade = case_when(Average_mark >= 70 ~ "First Class",
Average_mark >= 60 ~ "Upper Second Class Honour",
Average_mark >= 50 ~ "Lower Second Class Honour",
Average_mark >= 40 ~ "Third class/pass",
Average_mark < 40 ~ "Fail")) %>%
group_by(Ethnicity, Grade) %>%
summarise(counts = n()) %>%
mutate(freq = (counts / sum(counts))) %>%
mutate(percentage = scales::label_percent()(counts / sum(counts))) %>%
ggplot(aes(x=Ethnicity, y = freq, fill =Grade)) +
geom_bar(stat = 'identity') +
coord_flip() +
scale_fill_brewer(palette = "Pastel1") +
geom_text_repel(aes(label = percentage), position = position_stack(vjust = 0.5), size = 2.5) +
ggtitle("Distribution of degree outcomes by ethnic group (2017-2021)") +
theme(plot.title = element_text(hjust = 0.5)) +
labs(y = "Frequency", fill = "Degree Classification") #+
`summarise()` has grouped output by 'Ethnicity'. You can override using the `.groups` argument.
#scale_fill_discrete(breaks=c('First Class', 'Upper Second Class Honour', 'Lower Second Class Honour', 'Third class/pass', 'Fail'))
I have began by visualising the average mark for each assessment type (2017-2021)
#Visualise effect of assessment type on average mark between 2017-2021
full_data %>%
select(ASSESSMENT_COMPONENT_TYPE_DESC, MARK) %>%
group_by(ASSESSMENT_COMPONENT_TYPE_DESC) %>%
tally(mean(MARK)) %>%
rename(Average_mark = n) %>%
#spread(key = Ethnicity, value = Average_mark) %>%
ggplot(aes(x=ASSESSMENT_COMPONENT_TYPE_DESC, y = Average_mark, fill = ASSESSMENT_COMPONENT_TYPE_DESC)) +
geom_bar(stat = 'identity', position = 'dodge') +
scale_fill_brewer(palette = "Pastel2") +
geom_text(aes(label = round(Average_mark, 2)), position=position_dodge(width=0.9), vjust=-0.25) +
ggtitle("Average Mark for Coursework vs Exam (2017-2021)") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab('Assessment Type') +
ylab('Average Mark') +
guides(fill = guide_legend(title = "Assessment type"))
This bar graph above is showing that students overall perform much better in coursework than exams…but does ethnicity have an effect on this?? Will look at this later
Next i wanted to evaluate the average yearly mark dependent on assessment type. This shows every year the average grade is much bigger for cw type assessments vs exams. Essentially students perform much better in coursework type assessments in comparison to exams.
full_data %>%
select(INSTANCE_ACADEMIC_YEAR, ASSESSMENT_COMPONENT_TYPE_DESC,MARK) %>%
group_by(INSTANCE_ACADEMIC_YEAR, ASSESSMENT_COMPONENT_TYPE_DESC) %>%
tally(mean(MARK)) %>%
rename(Average_mark = n) %>%
spread(key = ASSESSMENT_COMPONENT_TYPE_DESC, value = Average_mark) %>%
mutate(Gap = Coursework - Exam)
This line graph below shows many things: 1. students perform much better in cw vs exams 2. Average cw mark is reducing over the years 3. Average exam mark is increasing over the years, but goes down in 2021 4. assessment awarding gap is reducing.
#Visualisation of above
full_data %>%
select(INSTANCE_ACADEMIC_YEAR, ASSESSMENT_COMPONENT_TYPE_DESC,MARK) %>%
group_by(INSTANCE_ACADEMIC_YEAR, ASSESSMENT_COMPONENT_TYPE_DESC) %>%
tally(mean(MARK)) %>%
rename(Average_mark = n) %>%
#spread(key = ASSESSMENT_COMPONENT_TYPE_DESC, value = Average_mark) %>%
#mutate(Gap = Coursework - Exam) %>%
ggplot(aes(x = INSTANCE_ACADEMIC_YEAR, y = Average_mark, col = ASSESSMENT_COMPONENT_TYPE_DESC)) +
geom_line() +
geom_point() +
ggtitle("Average yearly Mark for Coursework vs Exam (2017-2021)") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab('Academic Year') +
ylab('Average Mark') +
labs(col='Assessment Type')
Analysis of effect of assessment type on achieving top degree(2017-2022)
What about effect of assessment type on top degrees? This bar chart shows the average number of students that did cw vs exam in all students that achieved a top degree.
#average percentage of students who did cw vs exam to get first
full_data %>%
group_by(ID2) %>%
mutate(Average_mark = mean(MARK)) %>%
filter(Average_mark >=60) %>%
select(ID2, ASSESSMENT_COMPONENT_TYPE_DESC, Average_mark) %>%
distinct() %>%
group_by(ASSESSMENT_COMPONENT_TYPE_DESC) %>%
summarise(counts = n()) %>%
mutate(percentage = scales::label_percent()(counts / sum(counts))) %>%
ggplot(aes(x=ASSESSMENT_COMPONENT_TYPE_DESC, y = counts, fill = ASSESSMENT_COMPONENT_TYPE_DESC)) +
geom_bar(stat = 'identity', position = 'dodge') +
scale_fill_brewer(palette = "Pastel2") +
geom_text(aes(label = percentage), position=position_dodge(width=0.9), vjust=-0.25) +
ggtitle("Percentage of students doing cw vs exam for top degree") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab('Assessment Type') +
ylab('Number of students') +
guides(fill = guide_legend(title = "Assessment type"))
This bar chart is showing yearly average percentage of students who did cw vs exam to get first This results are pretty consistent over the years. Has reduced slightly over the years. Gap in 21 looks pretty small, but number of students in this year was lower.
#percentage of students who did cw vs exam to get first (2017-2022)
full_data %>%
group_by(ID2) %>%
mutate(Average_mark = mean(MARK)) %>%
filter(Average_mark >=60) %>%
select(ID2, ASSESSMENT_COMPONENT_TYPE_DESC, Average_mark, INSTANCE_ACADEMIC_YEAR) %>%
distinct() %>%
group_by(INSTANCE_ACADEMIC_YEAR, ASSESSMENT_COMPONENT_TYPE_DESC) %>%
summarise(counts = n()) %>%
mutate(percentage = scales::label_percent()(counts / sum(counts))) %>%
ggplot(aes(x=INSTANCE_ACADEMIC_YEAR, y = counts, fill =ASSESSMENT_COMPONENT_TYPE_DESC)) +
geom_bar(stat = 'identity', position = 'dodge') +
scale_fill_brewer(palette = "Pastel1") +
geom_text(aes(label = percentage), position=position_dodge(width=0.9), vjust=-0.25) +
ggtitle("Percentage of students doing cw vs exam for top degree each year") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab('Academic Year') +
ylab('Number of students') +
guides(fill = guide_legend(title = "Assessment type"))
`summarise()` has grouped output by 'INSTANCE_ACADEMIC_YEAR'. You can override using the `.groups` argument.
#Checking there is a reasonable number of white/others in each assessment type.
full_data %>%
select(Ethnicity, ASSESSMENT_COMPONENT_TYPE_DESC) %>%
group_by(Ethnicity, ASSESSMENT_COMPONENT_TYPE_DESC) %>%
tally()
First i want to look at affect of assessment type on student performance(White vs Other) This table below shows the difference in avgerage cw and exam marks for white Vs Bame students
Average grade over past 5 years for assessment type (AT) and ethnicity
full_data %>%
select(Ethnicity, ASSESSMENT_COMPONENT_TYPE_DESC,MARK) %>%
group_by(Ethnicity, ASSESSMENT_COMPONENT_TYPE_DESC) %>%
tally(mean(MARK)) %>%
rename(Average_mark = n) %>%
spread(key = Ethnicity, value = Average_mark) %>%
mutate(gap = White - BAME)
The gap doesnt look huge, but it is bigger for coursework type assessments. However this is gap in average mark (student performance), not the number of students who achieved a top degree. Therefore, this is not relating to Bame awarding gap. Still interesting…
Visualisation of above table Looking at average grade over past 5 years for AT and ethnicity
full_data %>%
select(Ethnicity, ASSESSMENT_COMPONENT_TYPE_DESC,MARK) %>%
group_by(Ethnicity, ASSESSMENT_COMPONENT_TYPE_DESC) %>%
tally(mean(MARK)) %>%
rename(Average_mark = n) %>%
#spread(key = Ethnicity, value = Average_mark) %>%
ggplot(aes(x=ASSESSMENT_COMPONENT_TYPE_DESC, y = Average_mark, fill = Ethnicity)) +
geom_bar(stat = 'identity', position = 'dodge') +
scale_fill_brewer(palette = "Accent") +
geom_text(aes(label = round(Average_mark, 2)), position=position_dodge(width=0.9), vjust=-0.25) +
ggtitle("Average Mark By Assessment Type And Ethnicity") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab('Assessment Type')
ylab('Average Mark')
$y
[1] "Average Mark"
attr(,"class")
[1] "labels"
This graph above is useful as it shows difference in average mark for each assessment type by ethnicity. However this research question pertains to the affect of assessment type on BAME awarding gaps. Therefore, instead of average marks, we need to compare the number of white vs BAME students who achieved top degrees and see if gap is greater for cw type assessments or coursework.
full_data %>%
group_by(ID2) %>%
mutate(Average_mark = mean(MARK)) %>%
filter(Average_mark >=60) %>%
select(ID2, Ethnicity, ASSESSMENT_COMPONENT_TYPE_DESC, Average_mark) %>%
distinct() %>%
select(Ethnicity, ASSESSMENT_COMPONENT_TYPE_DESC) %>%
group_by(Ethnicity,ASSESSMENT_COMPONENT_TYPE_DESC) %>%
summarise(counts = n()) %>%
mutate(percentage = scales::label_percent()(counts / sum(counts))) %>%
#select(Ethnicity, ASSESSMENT_COMPONENT_TYPE_DESC, percentage) %>%
ggplot(aes(x=ASSESSMENT_COMPONENT_TYPE_DESC, y = percentage, fill = Ethnicity)) +
geom_bar(stat = 'identity', position = 'dodge') +
scale_fill_brewer(palette = "Set1") +
geom_text(aes(label = percentage), position=position_dodge(width=0.9), vjust=-0.25) +
ggtitle("Effect of assessment type on awarding gap") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab('Assessment Type') +
ylab('Percentage of students')
Adding missing grouping variables: `ID2`
`summarise()` has grouped output by 'Ethnicity'. You can override using the `.groups` argument.
The results here are interesting. There is a 4% awarding gap for both ethnicities in each assessment type- this suggests exam do not lead to a wider awarding gap for the years 2017-2021.White students still performing better for both assessment types, but both ethnicities performing similar for assessment type.
#interaction plot of full dataset.
interaction.plot(x.factor = full_data$Ethnicity, #x-axis variable
trace.factor = full_data$ASSESSMENT_COMPONENT_TYPE_DESC, #variable for lines
response = full_data$MARK, #y-axis variable
fun = median, #metric to plot
ylab = "Average Mark",
xlab = "Ethnicity",
col = c("red", "blue"),
lty = 1, #line type
lwd = 2, #line width
trace.label = "Assessment type")
#Question 1:
# 1. What affect does the assessment type have on student performance?
#fitting model
mmodel1 <- lmer(MARK ~ ASSESSMENT_COMPONENT_TYPE_DESC * Ethnicity + PROGRAMME_PART + (1|ID2) + (1|MODULE_CODE), data = full_data)
summary(mmodel1)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: MARK ~ ASSESSMENT_COMPONENT_TYPE_DESC * Ethnicity + PROGRAMME_PART +
(1 | ID2) + (1 | MODULE_CODE)
Data: full_data
REML criterion at convergence: 6688052
Scaled residuals:
Min 1Q Median 3Q Max
-7.4683 -0.5331 0.0220 0.5787 4.7933
Random effects:
Groups Name Variance Std.Dev.
ID2 (Intercept) 38.83 6.232
MODULE_CODE (Intercept) 25.44 5.043
Residual 153.51 12.390
Number of obs: 840368, groups: ID2, 31540; MODULE_CODE, 6945
Fixed effects:
Estimate Std. Error df t value
(Intercept) 6.537e+01 1.521e-01 1.530e+04 429.825
ASSESSMENT_COMPONENT_TYPE_DESCExam -8.123e+00 5.911e-02 7.942e+05 -137.425
EthnicityWhite 3.333e+00 9.279e-02 3.443e+04 35.920
PROGRAMME_PARTB 5.245e-02 1.613e-01 1.550e+04 0.325
PROGRAMME_PARTC 7.652e-01 1.897e-01 1.295e+04 4.033
PROGRAMME_PARTD -1.309e+00 3.776e-01 3.151e+04 -3.467
PROGRAMME_PARTF 5.816e+00 5.340e-01 6.400e+03 10.892
PROGRAMME_PARTT -3.757e-01 1.903e-01 1.176e+04 -1.974
ASSESSMENT_COMPONENT_TYPE_DESCExam:EthnicityWhite 3.466e-01 6.729e-02 8.341e+05 5.151
Pr(>|t|)
(Intercept) < 2e-16 ***
ASSESSMENT_COMPONENT_TYPE_DESCExam < 2e-16 ***
EthnicityWhite < 2e-16 ***
PROGRAMME_PARTB 0.744988
PROGRAMME_PARTC 5.53e-05 ***
PROGRAMME_PARTD 0.000527 ***
PROGRAMME_PARTF < 2e-16 ***
PROGRAMME_PARTT 0.048427 *
ASSESSMENT_COMPONENT_TYPE_DESCExam:EthnicityWhite 2.59e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) ASSESSMENT_COMPONENT_TYPE_DESCEx EthncW
ASSESSMENT_COMPONENT_TYPE_DESCEx -0.130
EthnictyWht -0.414 0.155
PROGRAMME_PARTB -0.548 -0.001 0.003
PROGRAMME_PARTC -0.525 0.000 0.006
PROGRAMME_PARTD -0.281 0.004 0.011
PROGRAMME_PARTF -0.216 0.003 0.008
PROGRAMME_PARTT -0.716 0.037 0.155
ASSESSMENT_COMPONENT_TYPE_DESCE: 0.083 -0.709 -0.214
PROGRAMME_PARTB PROGRAMME_PARTC PROGRAMME_PARTD
ASSESSMENT_COMPONENT_TYPE_DESCEx
EthnictyWht
PROGRAMME_PARTB
PROGRAMME_PARTC 0.510
PROGRAMME_PARTD 0.245 0.328
PROGRAMME_PARTF 0.153 0.145 0.075
PROGRAMME_PARTT 0.440 0.427 0.255
ASSESSMENT_COMPONENT_TYPE_DESCE: -0.003 0.000 0.000
PROGRAMME_PARTF PROGRAMME_PARTT
ASSESSMENT_COMPONENT_TYPE_DESCEx
EthnictyWht
PROGRAMME_PARTB
PROGRAMME_PARTC
PROGRAMME_PARTD
PROGRAMME_PARTF
PROGRAMME_PARTT 0.170
ASSESSMENT_COMPONENT_TYPE_DESCE: -0.001 -0.010
optimizer (nloptwrap) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.00454631 (tol = 0.002, component 1)
#interaction plot of model...
plot_model(mmodel1, type = "pred", terms = c("Ethnicity", "ASSESSMENT_COMPONENT_TYPE_DESC"))
anova(mmodel1)
Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
ASSESSMENT_COMPONENT_TYPE_DESC 5366486 5366486 1 717508 34959.266 < 2.2e-16 ***
Ethnicity 224560 224560 1 32669 1462.867 < 2.2e-16 ***
PROGRAMME_PART 27274 5455 5 12989 35.534 < 2.2e-16 ***
ASSESSMENT_COMPONENT_TYPE_DESC:Ethnicity 4073 4073 1 834072 26.532 2.593e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#shows significance of all these variables on mark??
Results from the anova test show each of these variables have significant effect on mark
#effect size...look at video for how to explain this in report.
effectsize::eta_squared(mmodel1, partial = TRUE)
# Effect Size for ANOVA (Type III)
Parameter | Eta2 (partial) | 95% CI
------------------------------------------------------------------------
ASSESSMENT_COMPONENT_TYPE_DESC | 0.05 | [0.05, 1.00]
Ethnicity | 0.04 | [0.04, 1.00]
PROGRAMME_PART | 0.01 | [0.01, 1.00]
ASSESSMENT_COMPONENT_TYPE_DESC:Ethnicity | 3.18e-05 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
#emmeans is meant to be looked at for main effect.
#in this case they were all significant so will do for each??
emmeans(mmodel1, list(pairwise~ASSESSMENT_COMPONENT_TYPE_DESC), adjust="tukey")
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'pbkrtest.limit = 840368' (or larger)
[or, globally, 'set emm_options(pbkrtest.limit = 840368)' or larger];
but be warned that this may result in large computation time and memory use.
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'lmerTest.limit = 840368' (or larger)
[or, globally, 'set emm_options(lmerTest.limit = 840368)' or larger];
but be warned that this may result in large computation time and memory use.
NOTE: Results may be misleading due to involvement in interactions
$`emmeans of ASSESSMENT_COMPONENT_TYPE_DESC`
ASSESSMENT_COMPONENT_TYPE_DESC emmean SE df asymp.LCL asymp.UCL
Coursework 67.86 0.1263 Inf 67.62 68.11
Exam 59.91 0.1290 Inf 59.66 60.17
Results are averaged over the levels of: Ethnicity, PROGRAMME_PART
Degrees-of-freedom method: asymptotic
Confidence level used: 0.95
$`pairwise differences of ASSESSMENT_COMPONENT_TYPE_DESC`
1 estimate SE df z.ratio p.value
Coursework - Exam 7.95 0.0425 Inf 186.974 <.0001
Results are averaged over the levels of: Ethnicity, PROGRAMME_PART
Degrees-of-freedom method: asymptotic
#here comparison shows that there is a sig diff in mark when assessment type is cw compared to exam (p.value = <0.0001)
VarCorr(mmodel1)
Groups Name Std.Dev.
ID2 (Intercept) 6.2315
MODULE_CODE (Intercept) 5.0433
Residual 12.3898
totSD <- sqrt( 6.2315 + 5.0433 + 12.3898)
emmV <- emmeans(mmodel1, ~ ASSESSMENT_COMPONENT_TYPE_DESC)
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'pbkrtest.limit = 840368' (or larger)
[or, globally, 'set emm_options(pbkrtest.limit = 840368)' or larger];
but be warned that this may result in large computation time and memory use.
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'lmerTest.limit = 840368' (or larger)
[or, globally, 'set emm_options(lmerTest.limit = 840368)' or larger];
but be warned that this may result in large computation time and memory use.
NOTE: Results may be misleading due to involvement in interactions
eff_size(emmV, sigma = totSD, edf = 5)
contrast effect.size SE df asymp.LCL asymp.UCL
Coursework - Exam 1.63 0.517 Inf 0.621 2.65
Results are averaged over the levels of: Ethnicity, PROGRAMME_PART
sigma used for effect sizes: 4.865
Degrees-of-freedom method: inherited from asymptotic when re-gridding
Confidence level used: 0.95
emmeans(mmodel1, list(pairwise~Ethnicity), adjust="tukey")
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'pbkrtest.limit = 840368' (or larger)
[or, globally, 'set emm_options(pbkrtest.limit = 840368)' or larger];
but be warned that this may result in large computation time and memory use.
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'lmerTest.limit = 840368' (or larger)
[or, globally, 'set emm_options(lmerTest.limit = 840368)' or larger];
but be warned that this may result in large computation time and memory use.
NOTE: Results may be misleading due to involvement in interactions
$`emmeans of Ethnicity`
Ethnicity emmean SE df asymp.LCL asymp.UCL
BAME 62.14 0.1373 Inf 61.87 62.41
White 65.64 0.1305 Inf 65.39 65.90
Results are averaged over the levels of: ASSESSMENT_COMPONENT_TYPE_DESC, PROGRAMME_PART
Degrees-of-freedom method: asymptotic
Confidence level used: 0.95
$`pairwise differences of Ethnicity`
1 estimate SE df z.ratio p.value
BAME - White -3.51 0.0917 Inf -38.247 <.0001
Results are averaged over the levels of: ASSESSMENT_COMPONENT_TYPE_DESC, PROGRAMME_PART
Degrees-of-freedom method: asymptotic
#same interpretation as above but for ethnicity??? but ask safa why is other-white not white-other.
emmV2 <- emmeans(mmodel1, ~ Ethnicity)
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'pbkrtest.limit = 840368' (or larger)
[or, globally, 'set emm_options(pbkrtest.limit = 840368)' or larger];
but be warned that this may result in large computation time and memory use.
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'lmerTest.limit = 840368' (or larger)
[or, globally, 'set emm_options(lmerTest.limit = 840368)' or larger];
but be warned that this may result in large computation time and memory use.
NOTE: Results may be misleading due to involvement in interactions
eff_size(emmV2, sigma = totSD, edf = 5)
contrast effect.size SE df asymp.LCL asymp.UCL
BAME - White -0.721 0.229 Inf -1.17 -0.273
Results are averaged over the levels of: ASSESSMENT_COMPONENT_TYPE_DESC, PROGRAMME_PART
sigma used for effect sizes: 4.865
Degrees-of-freedom method: inherited from asymptotic when re-gridding
Confidence level used: 0.95
emmeans(mmodel1, list(pairwise~PROGRAMME_PART), adjust="tukey")
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'pbkrtest.limit = 840368' (or larger)
[or, globally, 'set emm_options(pbkrtest.limit = 840368)' or larger];
but be warned that this may result in large computation time and memory use.
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'lmerTest.limit = 840368' (or larger)
[or, globally, 'set emm_options(lmerTest.limit = 840368)' or larger];
but be warned that this may result in large computation time and memory use.
$`emmeans of PROGRAMME_PART`
PROGRAMME_PART emmean SE df asymp.LCL asymp.UCL
A 63.06 0.1387 Inf 62.79 63.34
B 63.12 0.1356 Inf 62.85 63.38
C 63.83 0.1582 Inf 63.52 64.14
D 61.76 0.3606 Inf 61.05 62.46
F 68.88 0.5194 Inf 67.86 69.90
T 62.69 0.1308 Inf 62.43 62.95
Results are averaged over the levels of: ASSESSMENT_COMPONENT_TYPE_DESC, Ethnicity
Degrees-of-freedom method: asymptotic
Confidence level used: 0.95
$`pairwise differences of PROGRAMME_PART`
1 estimate SE df z.ratio p.value
A - B -0.0525 0.161 Inf -0.325 0.9995
A - C -0.7652 0.190 Inf -4.033 0.0008
A - D 1.3093 0.378 Inf 3.467 0.0070
A - F -5.8159 0.534 Inf -10.892 <.0001
A - T 0.3757 0.190 Inf 1.974 0.3574
B - C -0.7127 0.175 Inf -4.062 0.0007
B - D 1.3617 0.372 Inf 3.656 0.0035
B - F -5.7634 0.534 Inf -10.800 <.0001
B - T 0.4281 0.188 Inf 2.282 0.2014
C - D 2.0744 0.363 Inf 5.717 <.0001
C - F -5.0507 0.540 Inf -9.350 <.0001
C - T 1.1408 0.203 Inf 5.609 <.0001
D - F -7.1252 0.630 Inf -11.304 <.0001
D - T -0.9336 0.377 Inf -2.476 0.1312
F - T 6.1916 0.536 Inf 11.561 <.0001
Results are averaged over the levels of: ASSESSMENT_COMPONENT_TYPE_DESC, Ethnicity
Degrees-of-freedom method: asymptotic
P value adjustment: tukey method for comparing a family of 6 estimates
#interesting as only some parts show sig difference on mark... some parts e.g. B,T, show no sig diff on mark?? wonder why?? is this correct interpretation??
emmeans(mmodel1, pairwise ~ Ethnicity | ASSESSMENT_COMPONENT_TYPE_DESC)
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'pbkrtest.limit = 840368' (or larger)
[or, globally, 'set emm_options(pbkrtest.limit = 840368)' or larger];
but be warned that this may result in large computation time and memory use.
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'lmerTest.limit = 840368' (or larger)
[or, globally, 'set emm_options(lmerTest.limit = 840368)' or larger];
but be warned that this may result in large computation time and memory use.
$emmeans
ASSESSMENT_COMPONENT_TYPE_DESC = Coursework:
Ethnicity emmean SE df asymp.LCL asymp.UCL
BAME 66.20 0.1380 Inf 65.93 66.47
White 69.53 0.1311 Inf 69.27 69.79
ASSESSMENT_COMPONENT_TYPE_DESC = Exam:
Ethnicity emmean SE df asymp.LCL asymp.UCL
BAME 58.07 0.1428 Inf 57.79 58.35
White 61.75 0.1345 Inf 61.49 62.02
Results are averaged over the levels of: PROGRAMME_PART
Degrees-of-freedom method: asymptotic
Confidence level used: 0.95
$contrasts
ASSESSMENT_COMPONENT_TYPE_DESC = Coursework:
contrast estimate SE df z.ratio p.value
BAME - White -3.33 0.0928 Inf -35.920 <.0001
ASSESSMENT_COMPONENT_TYPE_DESC = Exam:
contrast estimate SE df z.ratio p.value
BAME - White -3.68 0.1023 Inf -35.974 <.0001
Results are averaged over the levels of: PROGRAMME_PART
Degrees-of-freedom method: asymptotic
emmV3 <- emmeans(mmodel1, ~ Ethnicity | ASSESSMENT_COMPONENT_TYPE_DESC)
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'pbkrtest.limit = 840368' (or larger)
[or, globally, 'set emm_options(pbkrtest.limit = 840368)' or larger];
but be warned that this may result in large computation time and memory use.
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'lmerTest.limit = 840368' (or larger)
[or, globally, 'set emm_options(lmerTest.limit = 840368)' or larger];
but be warned that this may result in large computation time and memory use.
eff_size(emmV3, sigma = totSD, edf = 5)
ASSESSMENT_COMPONENT_TYPE_DESC = Coursework:
contrast effect.size SE df asymp.LCL asymp.UCL
BAME - White -0.685 0.218 Inf -1.11 -0.259
ASSESSMENT_COMPONENT_TYPE_DESC = Exam:
contrast effect.size SE df asymp.LCL asymp.UCL
BAME - White -0.756 0.240 Inf -1.23 -0.286
Results are averaged over the levels of: PROGRAMME_PART
sigma used for effect sizes: 4.865
Degrees-of-freedom method: inherited from asymptotic when re-gridding
Confidence level used: 0.95
Usually, there are 3 model assumptions which we need to check; linearity, homoscedasticity and normality of residuals. We do not need to look at linearity as all variables in this model are categorical.
#code adapted from https://ademos.people.uic.edu/Chapter18.html
plot(mmodel1)
boxplot(mmodel1@resp$wtres~full_data$ASSESSMENT_COMPONENT_TYPE_DESC)
boxplot(mmodel1@resp$wtres~full_data$Ethnicity)
boxplot(mmodel1@resp$wtres~full_data$PROGRAMME_PART)
Another way to check for homoscedasticity - (sd) residuals
#ASSESSMENT TYPE
mmodel1res<-mmodel1@resp$wtres
#assessment type residuals in df
mmodel1df<-data.frame(full_data$ASSESSMENT_COMPONENT_TYPE_DESC,mmodel1res)
names(mmodel1df)<-c("ASSESSMENT_COMPONENT_TYPE_DESC","residual")
mmodel1df %>% group_by(ASSESSMENT_COMPONENT_TYPE_DESC)%>% summarise(sd(residual))
Since the highest sd is less than twice the lowest, that is fine (according to general rule of thumb for homoscedasticity across groups of a factor)
Now we will check for the remaining variables…
#Ethnicity
mmodel1df2<-data.frame(full_data$Ethnicity,mmodel1res)
names(mmodel1df2)<-c("Ethnicity","residual")
mmodel1df2 %>% group_by(Ethnicity)%>% summarise(sd(residual))
#assumption satisfied
#programme part
mmodel1df3<-data.frame(full_data$PROGRAMME_PART,mmodel1res)
names(mmodel1df3)<-c("Programme_Part","residual")
mmodel1df3 %>% group_by(Programme_Part)%>% summarise(sd(residual))
#assumption met for this categorical variable as the largest residual is 15.3 and is less than 2* the smallest residual of D (9.20)
qqnorm(mmodel1res)
#Interval estimates for mixed models
#confint(mmodel1)
The above model is assessing effect of variables on student performance by looking at average marks. We want to look at awarding gaps so will categorise the marks into top degree and below and run a binomial model.
#creating binomial response variable
full_data %>%
group_by(ID2) %>%
mutate(Average_mark = mean(MARK)) %>%
mutate(Grade = case_when(Average_mark >= 60 ~ 1,
Average_mark < 60 ~ 0)) -> full_data3
#1. What affect does the assessment type have on the BAME awarding gap?
linmod<-glmer(as.factor(Grade) ~ ASSESSMENT_COMPONENT_TYPE_DESC * Ethnicity + PROGRAMME_PART + (1|ID2) + (1|MODULE_CODE), family="binomial",data = full_data3, nAGQ=0)
summary(linmod)
Generalized linear mixed model fit by maximum likelihood (Adaptive Gauss-Hermite Quadrature,
nAGQ = 0) [glmerMod]
Family: binomial ( logit )
Formula: as.factor(Grade) ~ ASSESSMENT_COMPONENT_TYPE_DESC * Ethnicity +
PROGRAMME_PART + (1 | ID2) + (1 | MODULE_CODE)
Data: full_data3
AIC BIC logLik deviance df.resid
52332.8 52460.9 -26155.4 52310.8 840357
Scaled residuals:
Min 1Q Median 3Q Max
-0.18149 0.01053 0.01402 0.01861 0.07701
Random effects:
Groups Name Variance Std.Dev.
ID2 (Intercept) 3.070e+02 1.752e+01
MODULE_CODE (Intercept) 6.814e-12 2.610e-06
Number of obs: 840368, groups: ID2, 31540; MODULE_CODE, 6945
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.5760 0.2312 15.466 < 2e-16 ***
ASSESSMENT_COMPONENT_TYPE_DESCExam -0.1197 0.1662 -0.720 0.47152
EthnicityWhite 2.3745 0.2582 9.196 < 2e-16 ***
PROGRAMME_PARTB 0.4248 0.1625 2.614 0.00895 **
PROGRAMME_PARTC 0.5394 0.2511 2.148 0.03170 *
PROGRAMME_PARTD 1.9214 1.8806 1.022 0.30694
PROGRAMME_PARTF -0.5267 0.5435 -0.969 0.33257
PROGRAMME_PARTT -0.1385 0.2593 -0.534 0.59341
ASSESSMENT_COMPONENT_TYPE_DESCExam:EthnicityWhite -0.0483 0.2357 -0.205 0.83768
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) ASSESSMENT_COMPONENT_TYPE_DESCEx EthncW
ASSESSMENT_COMPONENT_TYPE_DESCEx -0.237
EthnictyWht -0.728 0.196
PROGRAMME_PARTB -0.185 -0.019 0.006
PROGRAMME_PARTC -0.134 -0.009 0.007
PROGRAMME_PARTD -0.020 0.003 0.002
PROGRAMME_PARTF -0.131 0.008 0.030
PROGRAMME_PARTT -0.687 0.095 0.392
ASSESSMENT_COMPONENT_TYPE_DESCE: 0.144 -0.701 -0.272
PROGRAMME_PARTB PROGRAMME_PARTC PROGRAMME_PARTD
ASSESSMENT_COMPONENT_TYPE_DESCEx
EthnictyWht
PROGRAMME_PARTB
PROGRAMME_PARTC 0.282
PROGRAMME_PARTD 0.037 0.038
PROGRAMME_PARTF 0.080 0.052 0.008
PROGRAMME_PARTT 0.168 0.122 0.017
ASSESSMENT_COMPONENT_TYPE_DESCE: -0.001 0.000 -0.002
PROGRAMME_PARTF PROGRAMME_PARTT
ASSESSMENT_COMPONENT_TYPE_DESCEx
EthnictyWht
PROGRAMME_PARTB
PROGRAMME_PARTC
PROGRAMME_PARTD
PROGRAMME_PARTF
PROGRAMME_PARTT 0.107
ASSESSMENT_COMPONENT_TYPE_DESCE: 0.001 -0.030
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
se <- sqrt(diag(vcov(linmod)))
# table of estimates with 95% CI
(tab <- cbind(Est = fixef(linmod), LL = fixef(linmod) - 1.96 * se, UL = fixef(linmod) + 1.96 *
se))
Est LL UL
(Intercept) 3.57595111 3.12276002 4.0291422
ASSESSMENT_COMPONENT_TYPE_DESCExam -0.11965227 -0.44536946 0.2060649
EthnicityWhite 2.37450940 1.86844174 2.8805771
PROGRAMME_PARTB 0.42478261 0.10625209 0.7433131
PROGRAMME_PARTC 0.53936330 0.04726069 1.0314659
PROGRAMME_PARTD 1.92137632 -1.76467657 5.6074292
PROGRAMME_PARTF -0.52666637 -1.59200991 0.5386772
PROGRAMME_PARTT -0.13846290 -0.64677264 0.3698468
ASSESSMENT_COMPONENT_TYPE_DESCExam:EthnicityWhite -0.04829525 -0.51035739 0.4137669
exp(tab)
Est LL UL
(Intercept) 35.7285865 22.7089705 56.212671
ASSESSMENT_COMPONENT_TYPE_DESCExam 0.8872289 0.6405876 1.228833
EthnicityWhite 10.7457400 6.4781938 17.824556
PROGRAMME_PARTB 1.5292579 1.1121022 2.102891
PROGRAMME_PARTC 1.7149146 1.0483953 2.805175
PROGRAMME_PARTD 6.8303527 0.1712422 272.442940
PROGRAMME_PARTF 0.5905704 0.2035162 1.713738
PROGRAMME_PARTT 0.8706956 0.5237333 1.447513
ASSESSMENT_COMPONENT_TYPE_DESCExam:EthnicityWhite 0.9528524 0.6002810 1.512504
#evaluating model significance/quality
#Analysis of Deviance
#install.packages("car")
#library("car")
Anova(linmod)
Analysis of Deviance Table (Type II Wald chisquare tests)
Response: as.factor(Grade)
Chisq Df Pr(>Chisq)
ASSESSMENT_COMPONENT_TYPE_DESC 1.4648 1 0.22616
Ethnicity 90.2155 1 < 2e-16 ***
PROGRAMME_PART 12.3015 5 0.03088 *
ASSESSMENT_COMPONENT_TYPE_DESC:Ethnicity 0.0420 1 0.83768
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
plot_model(linmod, type = "pred", terms = c("Ethnicity", "ASSESSMENT_COMPONENT_TYPE_DESC"))
emmeans(linmod, ~ ASSESSMENT_COMPONENT_TYPE_DESC, type = "response")
NOTE: Results may be misleading due to involvement in interactions
ASSESSMENT_COMPONENT_TYPE_DESC response SE df asymp.LCL asymp.UCL
Coursework 0.9941 0.002031 Inf 0.9885 0.9970
Exam 0.9932 0.002389 Inf 0.9865 0.9966
Results are averaged over the levels of: Ethnicity, PROGRAMME_PART
Unknown transformation "as.factor": no transformation done
Confidence level used: 0.95
emmeans(linmod, ~ Ethnicity, type = "response")
NOTE: Results may be misleading due to involvement in interactions
Ethnicity response SE df asymp.LCL asymp.UCL
BAME 0.9799 0.0074121 Inf 0.9589 0.9903
White 0.9980 0.0007064 Inf 0.9960 0.9990
Results are averaged over the levels of: ASSESSMENT_COMPONENT_TYPE_DESC, PROGRAMME_PART
Unknown transformation "as.factor": no transformation done
Confidence level used: 0.95
emmeans(linmod, ~ Ethnicity | ASSESSMENT_COMPONENT_TYPE_DESC, type = "response")
ASSESSMENT_COMPONENT_TYPE_DESC = Coursework:
Ethnicity response SE df asymp.LCL asymp.UCL
BAME 0.9810 0.0070326 Inf 0.9610 0.9909
White 0.9982 0.0006544 Inf 0.9963 0.9991
ASSESSMENT_COMPONENT_TYPE_DESC = Exam:
Ethnicity response SE df asymp.LCL asymp.UCL
BAME 0.9787 0.0081866 Inf 0.9551 0.9900
White 0.9979 0.0008030 Inf 0.9955 0.9990
Results are averaged over the levels of: PROGRAMME_PART
Unknown transformation "as.factor": no transformation done
Confidence level used: 0.95
This model does not require checking of assumptions
Data wrangling: This research question required different categorisisation of ethnicities.
#categorising ethnicities into Asian, Black, Chinese, Mixed & White
clean_ethnicity_data2[clean_ethnicity_data2 == "White - British"] <- "White"
clean_ethnicity_data2[clean_ethnicity_data2 == "White - Irish"] <- "White"
clean_ethnicity_data2[clean_ethnicity_data2 == "White - Scottish"] <- "White"
clean_ethnicity_data2[clean_ethnicity_data2 == "Other White Background"] <- "White"
clean_ethnicity_data2[clean_ethnicity_data2 == "Asian Or Asian British - Bangladeshi"] <- "Asian"
clean_ethnicity_data2[clean_ethnicity_data2 == "Asian Or Asian British - Indian"] <- "Asian"
clean_ethnicity_data2[clean_ethnicity_data2 == "Asian Or Asian British - Pakistani"] <- "Asian"
clean_ethnicity_data2[clean_ethnicity_data2 == "Black Or Black British - African"] <- "Black"
clean_ethnicity_data2[clean_ethnicity_data2 == "Black Or Black British - Caribbean"] <- "Black"
clean_ethnicity_data2[clean_ethnicity_data2 == "Chinese"] <- "Chinese"
clean_ethnicity_data2[clean_ethnicity_data2 == "Mixed - White And Asian"] <- "Mixed"
clean_ethnicity_data2[clean_ethnicity_data2 == "Mixed - White And Black African"] <- "Mixed"
clean_ethnicity_data2[clean_ethnicity_data2 == "Mixed - White And Black Caribbean"] <- "Mixed"
clean_ethnicity_data2[clean_ethnicity_data2 == "Other Asian Background"] <- "Asian"
clean_ethnicity_data2[clean_ethnicity_data2 == "Other Black Background"] <- "Black"
clean_ethnicity_data2[clean_ethnicity_data2 == "Other Mixed Background"] <- "Mixed"
#remove not known,no information & info refused- as could be any ethnicity.
#also removed gypsy as only 1 student- not good for representation
#also removed other ethnic background- as not sure what category this relates to
clean_ethnicity_data2<-clean_ethnicity_data2[!(clean_ethnicity_data2$Ethnicity=="Information Refused" | clean_ethnicity_data2$Ethnicity=="Information Refused - OBSOLETE VALUE" | clean_ethnicity_data2$Ethnicity=="No Information" | clean_ethnicity_data2$Ethnicity=="Not Known" | clean_ethnicity_data2$Ethnicity=="Gypsy or Traveller" | clean_ethnicity_data2$Ethnicity=="Other Ethnic Background" | clean_ethnicity_data2$Ethnicity=="Arab") ,]
clean_ethnicity_data2 %>%
mutate_at(vars(Gender, Ethnicity), list(factor)) -> clean_ethnicity_data2
summary(clean_ethnicity_data2)
ID2 Gender Ethnicity
Min. : 1 Man :24318 Asian : 5408
1st Qu.:11056 Woman:17820 Black : 2814
Median :22126 Chinese: 6688
Mean :22104 Mixed : 1865
3rd Qu.:33135 White :25363
Max. :44191
#Ethnicity
#Percentage of students belonging to each ethnicity after categorisation
clean_ethnicity_data2 %>%
group_by(Ethnicity) %>%
summarise(counts = n()) %>%
mutate(freq = (counts / sum(counts))) %>%
mutate(percentage = scales::label_percent()(counts / sum(counts))) %>%
ggplot( ., aes(x="", y=freq, fill=Ethnicity)) +
geom_bar(stat="identity", width=1, color="white") +
coord_polar("y", start=0) +
theme_void() +
geom_text(aes(label = percentage), color = "white", position = position_stack(vjust = 0.5), size = 2.5) +
scale_fill_brewer(palette = "Paired") +
ggtitle("Percentage of students in each ethnic category") +
theme(plot.title = element_text(hjust = 0.5))
Joining two datasets
#JOIN ETHNIC DATA AND MARK DATA
inner_join(clean_ethnicity_data2, clean_mark_data) -> full_data2
print(distinct(full_data2))
boxplot(MARK~Ethnicity, data = full_data2)
#Visualise effect of ethnicity on overall/mean mark between 2017-2021
full_data2 %>%
select(Ethnicity, MARK) %>%
group_by(Ethnicity) %>%
tally(mean(MARK)) %>%
rename(Average_mark = n) %>%
ggplot(aes(x=Ethnicity, y = Average_mark, fill =Ethnicity)) +
geom_bar(stat = 'identity', position = 'dodge') +
scale_fill_brewer(palette = "Paired") +
geom_text(aes(label = round(Average_mark, 2)), position=position_dodge(width=0.9), vjust=-0.25) +
ggtitle("Average Test Scores by Ethnicity between 2017-2021") +
theme(plot.title = element_text(hjust = 0.5)) +
ylab('Average Mark')
#looking at avg grade for each ethnicity type between 2017-2021
full_data2 %>%
select(INSTANCE_ACADEMIC_YEAR, Ethnicity,MARK) %>%
group_by(INSTANCE_ACADEMIC_YEAR, Ethnicity) %>%
tally(mean(MARK)) %>%
rename(Average_mark = n) %>%
#spread(key = Ethnicity, value = Average_mark) %>%
ggplot(aes(x = INSTANCE_ACADEMIC_YEAR, y = Average_mark, col = Ethnicity)) +
geom_line() +
geom_point() +
ggtitle("Average grade for White vs BAME (2017-2021)") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab('Academic Year') +
ylab('Average Mark') +
labs(col='Ethnicity')
#Visualising percentage of white vs BAME students achieving top degree (2017-2021)
full_data2 %>%
group_by(ID2) %>%
mutate(Average_mark = mean(MARK)) %>%
filter(Average_mark >=60) %>%
select(ID2, Ethnicity, Average_mark) %>%
distinct() %>%
group_by(Ethnicity) %>%
summarise(counts = n()) %>%
mutate(percentage = scales::label_percent()(counts / sum(counts))) %>%
ggplot(aes(x=Ethnicity, y = counts, fill =Ethnicity)) +
geom_bar(stat = 'identity', position = 'dodge') +
scale_fill_brewer(palette = "Pastel1") +
geom_text(aes(label = percentage), position = position_stack(vjust = 1.1), size = 3) +
ggtitle("Percentage of white vs BAME students achieving top degree (2017-2021)") +
theme(plot.title = element_text(hjust = 0.5)) +
ylab('Number of students')
#BAME AWARDING GAP
full_data2 %>%
group_by(ID2) %>%
mutate(Average_mark = mean(MARK)) %>%
mutate(Grade = case_when(Average_mark >= 70 ~ "First Class",
Average_mark >= 60 ~ "Upper Second Class Honour",
Average_mark >= 50 ~ "Lower Second Class Honour",
Average_mark >= 40 ~ "Third class/pass",
Average_mark < 40 ~ "Fail")) %>%
group_by(Ethnicity, Grade) %>%
summarise(counts = n()) %>%
mutate(freq = (counts / sum(counts))) %>%
mutate(percentage = scales::label_percent()(counts / sum(counts))) %>%
ggplot(aes(x=Ethnicity, y = freq, fill =Grade)) +
geom_bar(stat = 'identity') +
coord_flip() +
scale_fill_brewer(palette = "Pastel1") +
geom_text(aes(label = percentage), position = position_stack(vjust = 0.5), size = 2.5) +
ggtitle("Distribution of degree outcomes by ethnic group (2017-2021)") +
theme(plot.title = element_text(hjust = 0.5)) +
labs(y = "Frequency", fill = "Degree Classification")
`summarise()` has grouped output by 'Ethnicity'. You can override using the `.groups` argument.
full_data2 %>%
select(Ethnicity, ASSESSMENT_COMPONENT_TYPE_DESC,MARK) %>%
group_by(Ethnicity, ASSESSMENT_COMPONENT_TYPE_DESC) %>%
tally(mean(MARK)) %>%
rename(Average_mark = n) %>%
#spread(key = Ethnicity, value = Average_mark) %>%
ggplot(aes(x=ASSESSMENT_COMPONENT_TYPE_DESC, y = Average_mark, fill = Ethnicity)) +
geom_bar(stat = 'identity', position = 'dodge') +
scale_fill_brewer(palette = "Set1") +
geom_text(aes(label = round(Average_mark, 2)), position=position_dodge(width=0.9), vjust=-0.25) +
ggtitle("Average Mark By Assessment Type And Ethnicity") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab('Assessment Type') +
ylab('Average Mark')
Effect of assessment type on awarding gap for Non-white categorise.
full_data2 %>%
group_by(ID2) %>%
mutate(Average_mark = mean(MARK)) %>%
filter(Average_mark >=60) %>%
select(ID2, Ethnicity, ASSESSMENT_COMPONENT_TYPE_DESC, Average_mark) %>%
distinct() %>%
select(Ethnicity, ASSESSMENT_COMPONENT_TYPE_DESC) %>%
group_by(Ethnicity,ASSESSMENT_COMPONENT_TYPE_DESC) %>%
summarise(counts = n()) %>%
mutate(percentage = scales::label_percent()(counts / sum(counts))) %>%
#select(Ethnicity, ASSESSMENT_COMPONENT_TYPE_DESC, percentage) %>%
ggplot(aes(x=ASSESSMENT_COMPONENT_TYPE_DESC, y = percentage, fill = Ethnicity)) +
geom_bar(stat = 'identity', position = 'dodge') +
scale_fill_brewer(palette = "Set1") +
geom_text(aes(label = percentage), position=position_dodge(width=0.9), vjust=-0.25) +
ggtitle("Effect of assessment type on awarding gap") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab('Assessment Type') +
ylab('Percentage of students')
Adding missing grouping variables: `ID2`
`summarise()` has grouped output by 'Ethnicity'. You can override using the `.groups` argument.
interaction.plot(x.factor = full_data2$Ethnicity, #x-axis variable
trace.factor = full_data2$ASSESSMENT_COMPONENT_TYPE_DESC, #variable for lines
response = full_data2$MARK, #y-axis variable
fun = median, #metric to plot
ylab = "Mark",
xlab = "Ethnicity",
col = c("pink", "lightblue"),
lty = 1, #line type
lwd = 2, #line width
trace.label = "Assessment type")
#Question 2:
# 2. What affect does the assessment type have on the performance of Non-White students (Namely Asian, Black, Chinese, and mixed)
mmodel2<- lmer(MARK ~ ASSESSMENT_COMPONENT_TYPE_DESC * Ethnicity + PROGRAMME_PART + (1|ID2) + (1|MODULE_CODE), data = full_data2)
summary(mmodel2)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: MARK ~ ASSESSMENT_COMPONENT_TYPE_DESC * Ethnicity + PROGRAMME_PART +
(1 | ID2) + (1 | MODULE_CODE)
Data: full_data2
REML criterion at convergence: 6585014
Scaled residuals:
Min 1Q Median 3Q Max
-7.4776 -0.5326 0.0220 0.5783 4.8132
Random effects:
Groups Name Variance Std.Dev.
ID2 (Intercept) 37.93 6.158
MODULE_CODE (Intercept) 25.28 5.028
Residual 153.15 12.375
Number of obs: 827711, groups: ID2, 31062; MODULE_CODE, 6939
Fixed effects:
Estimate Std. Error df t value
(Intercept) 6.634e+01 1.747e-01 2.325e+04 379.791
ASSESSMENT_COMPONENT_TYPE_DESCExam -8.116e+00 8.779e-02 8.194e+05 -92.448
EthnicityBlack -1.485e+00 1.895e-01 3.176e+04 -7.835
EthnicityChinese -3.433e+00 1.631e-01 3.652e+04 -21.052
EthnicityMixed 9.747e-01 2.231e-01 3.176e+04 4.369
EthnicityWhite 2.288e+00 1.269e-01 3.332e+04 18.035
PROGRAMME_PARTB 4.924e-02 1.615e-01 1.522e+04 0.305
PROGRAMME_PARTC 7.788e-01 1.898e-01 1.277e+04 4.102
PROGRAMME_PARTD -1.202e+00 3.782e-01 3.099e+04 -3.179
PROGRAMME_PARTF 5.688e+00 5.341e-01 6.438e+03 10.649
PROGRAMME_PARTT 3.115e-01 1.935e-01 1.250e+04 1.609
ASSESSMENT_COMPONENT_TYPE_DESCExam:EthnicityBlack -1.281e+00 1.393e-01 8.183e+05 -9.197
ASSESSMENT_COMPONENT_TYPE_DESCExam:EthnicityChinese 1.423e+00 1.449e-01 8.248e+05 9.823
ASSESSMENT_COMPONENT_TYPE_DESCExam:EthnicityMixed -1.760e-01 1.620e-01 8.175e+05 -1.086
ASSESSMENT_COMPONENT_TYPE_DESCExam:EthnicityWhite 3.375e-01 9.320e-02 8.203e+05 3.622
Pr(>|t|)
(Intercept) < 2e-16 ***
ASSESSMENT_COMPONENT_TYPE_DESCExam < 2e-16 ***
EthnicityBlack 4.84e-15 ***
EthnicityChinese < 2e-16 ***
EthnicityMixed 1.25e-05 ***
EthnicityWhite < 2e-16 ***
PROGRAMME_PARTB 0.760414
PROGRAMME_PARTC 4.12e-05 ***
PROGRAMME_PARTD 0.001478 **
PROGRAMME_PARTF < 2e-16 ***
PROGRAMME_PARTT 0.107560
ASSESSMENT_COMPONENT_TYPE_DESCExam:EthnicityBlack < 2e-16 ***
ASSESSMENT_COMPONENT_TYPE_DESCExam:EthnicityChinese < 2e-16 ***
ASSESSMENT_COMPONENT_TYPE_DESCExam:EthnicityMixed 0.277308
ASSESSMENT_COMPONENT_TYPE_DESCExam:EthnicityWhite 0.000293 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation matrix not shown by default, as p = 15 > 12.
Use print(x, correlation=TRUE) or
vcov(x) if you need it
#same as above- but shows more values ...
plot_model(mmodel2, type = "pred", terms = c("Ethnicity", "ASSESSMENT_COMPONENT_TYPE_DESC"))
anova(mmodel2)
Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
ASSESSMENT_COMPONENT_TYPE_DESC 3439838 3439838 1 764529 22460.340 < 2.2e-16 ***
Ethnicity 277515 69379 4 32481 453.007 < 2.2e-16 ***
PROGRAMME_PART 23368 4674 5 13145 30.516 < 2.2e-16 ***
ASSESSMENT_COMPONENT_TYPE_DESC:Ethnicity 46918 11730 4 822301 76.588 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Similar to RQ1, all variables hold a significant effect on student performance
#effect size...look at video for how to explain this in report.
effectsize::eta_squared(mmodel2, partial = TRUE)
# Effect Size for ANOVA (Type III)
Parameter | Eta2 (partial) | 95% CI
------------------------------------------------------------------------
ASSESSMENT_COMPONENT_TYPE_DESC | 0.03 | [0.03, 1.00]
Ethnicity | 0.05 | [0.05, 1.00]
PROGRAMME_PART | 0.01 | [0.01, 1.00]
ASSESSMENT_COMPONENT_TYPE_DESC:Ethnicity | 3.72e-04 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
#emmeans is meant to be looked at for main effect.
#in this case they were all significant so will do for each??
emmeans(mmodel2, list(pairwise~ASSESSMENT_COMPONENT_TYPE_DESC), adjust="tukey")
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'pbkrtest.limit = 827711' (or larger)
[or, globally, 'set emm_options(pbkrtest.limit = 827711)' or larger];
but be warned that this may result in large computation time and memory use.
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'lmerTest.limit = 827711' (or larger)
[or, globally, 'set emm_options(lmerTest.limit = 827711)' or larger];
but be warned that this may result in large computation time and memory use.
NOTE: Results may be misleading due to involvement in interactions
$`emmeans of ASSESSMENT_COMPONENT_TYPE_DESC`
ASSESSMENT_COMPONENT_TYPE_DESC emmean SE df asymp.LCL asymp.UCL
Coursework 66.94 0.1335 Inf 66.68 67.20
Exam 58.89 0.1376 Inf 58.62 59.16
Results are averaged over the levels of: Ethnicity, PROGRAMME_PART
Degrees-of-freedom method: asymptotic
Confidence level used: 0.95
$`pairwise differences of ASSESSMENT_COMPONENT_TYPE_DESC`
1 estimate SE df z.ratio p.value
Coursework - Exam 8.06 0.0538 Inf 149.868 <.0001
Results are averaged over the levels of: Ethnicity, PROGRAMME_PART
Degrees-of-freedom method: asymptotic
emmeans(mmodel2, list(pairwise~Ethnicity), adjust="tukey")
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'pbkrtest.limit = 827711' (or larger)
[or, globally, 'set emm_options(pbkrtest.limit = 827711)' or larger];
but be warned that this may result in large computation time and memory use.
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'lmerTest.limit = 827711' (or larger)
[or, globally, 'set emm_options(lmerTest.limit = 827711)' or larger];
but be warned that this may result in large computation time and memory use.
NOTE: Results may be misleading due to involvement in interactions
$`emmeans of Ethnicity`
Ethnicity emmean SE df asymp.LCL asymp.UCL
Asian 63.22 0.1618 Inf 62.90 63.53
Black 61.09 0.1919 Inf 60.72 61.47
Chinese 60.49 0.1771 Inf 60.15 60.84
Mixed 64.10 0.2227 Inf 63.67 64.54
White 65.67 0.1303 Inf 65.42 65.93
Results are averaged over the levels of: ASSESSMENT_COMPONENT_TYPE_DESC, PROGRAMME_PART
Degrees-of-freedom method: asymptotic
Confidence level used: 0.95
$`pairwise differences of Ethnicity`
1 estimate SE df z.ratio p.value
Asian - Black 2.125 0.187 Inf 11.391 <.0001
Asian - Chinese 2.722 0.165 Inf 16.481 <.0001
Asian - Mixed -0.887 0.219 Inf -4.042 0.0005
Asian - White -2.457 0.125 Inf -19.692 <.0001
Black - Chinese 0.597 0.197 Inf 3.029 0.0207
Black - Mixed -3.012 0.243 Inf -12.409 <.0001
Black - White -4.582 0.163 Inf -28.163 <.0001
Chinese - Mixed -3.608 0.231 Inf -15.608 <.0001
Chinese - White -5.179 0.145 Inf -35.600 <.0001
Mixed - White -1.570 0.197 Inf -7.985 <.0001
Results are averaged over the levels of: ASSESSMENT_COMPONENT_TYPE_DESC, PROGRAMME_PART
Degrees-of-freedom method: asymptotic
P value adjustment: tukey method for comparing a family of 5 estimates
emmeans(mmodel2, list(pairwise~PROGRAMME_PART), adjust="tukey")
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'pbkrtest.limit = 827711' (or larger)
[or, globally, 'set emm_options(pbkrtest.limit = 827711)' or larger];
but be warned that this may result in large computation time and memory use.
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'lmerTest.limit = 827711' (or larger)
[or, globally, 'set emm_options(lmerTest.limit = 827711)' or larger];
but be warned that this may result in large computation time and memory use.
$`emmeans of PROGRAMME_PART`
PROGRAMME_PART emmean SE df asymp.LCL asymp.UCL
A 61.98 0.1461 Inf 61.69 62.26
B 62.03 0.1431 Inf 61.75 62.31
C 62.76 0.1646 Inf 62.43 63.08
D 60.78 0.3638 Inf 60.06 61.49
F 67.67 0.5214 Inf 66.64 68.69
T 62.29 0.1372 Inf 62.02 62.56
Results are averaged over the levels of: ASSESSMENT_COMPONENT_TYPE_DESC, Ethnicity
Degrees-of-freedom method: asymptotic
Confidence level used: 0.95
$`pairwise differences of PROGRAMME_PART`
1 estimate SE df z.ratio p.value
A - B -0.0492 0.161 Inf -0.305 0.9996
A - C -0.7788 0.190 Inf -4.102 0.0006
A - D 1.2023 0.378 Inf 3.179 0.0185
A - F -5.6881 0.534 Inf -10.649 <.0001
A - T -0.3115 0.194 Inf -1.609 0.5923
B - C -0.7296 0.176 Inf -4.151 0.0005
B - D 1.2515 0.373 Inf 3.355 0.0103
B - F -5.6388 0.534 Inf -10.563 <.0001
B - T -0.2623 0.191 Inf -1.374 0.7426
C - D 1.9811 0.363 Inf 5.450 <.0001
C - F -4.9093 0.540 Inf -9.085 <.0001
C - T 0.4673 0.206 Inf 2.266 0.2083
D - F -6.8904 0.631 Inf -10.923 <.0001
D - T -1.5138 0.379 Inf -3.993 0.0009
F - T 5.3766 0.537 Inf 10.012 <.0001
Results are averaged over the levels of: ASSESSMENT_COMPONENT_TYPE_DESC, Ethnicity
Degrees-of-freedom method: asymptotic
P value adjustment: tukey method for comparing a family of 6 estimates
#There was an interaction so need emmeans for this
emmeans(mmodel2, pairwise ~ Ethnicity | ASSESSMENT_COMPONENT_TYPE_DESC)
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'pbkrtest.limit = 827711' (or larger)
[or, globally, 'set emm_options(pbkrtest.limit = 827711)' or larger];
but be warned that this may result in large computation time and memory use.
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'lmerTest.limit = 827711' (or larger)
[or, globally, 'set emm_options(lmerTest.limit = 827711)' or larger];
but be warned that this may result in large computation time and memory use.
$emmeans
ASSESSMENT_COMPONENT_TYPE_DESC = Coursework:
Ethnicity emmean SE df asymp.LCL asymp.UCL
Asian 67.27 0.1633 Inf 66.95 67.59
Black 65.79 0.1939 Inf 65.41 66.17
Chinese 63.84 0.1750 Inf 63.50 64.18
Mixed 68.25 0.2255 Inf 67.81 68.69
White 69.56 0.1309 Inf 69.31 69.82
ASSESSMENT_COMPONENT_TYPE_DESC = Exam:
Ethnicity emmean SE df asymp.LCL asymp.UCL
Asian 59.16 0.1719 Inf 58.82 59.49
Black 56.39 0.2066 Inf 55.99 56.80
Chinese 57.15 0.1986 Inf 56.76 57.54
Mixed 59.96 0.2416 Inf 59.48 60.43
White 61.78 0.1343 Inf 61.52 62.05
Results are averaged over the levels of: PROGRAMME_PART
Degrees-of-freedom method: asymptotic
Confidence level used: 0.95
$contrasts
ASSESSMENT_COMPONENT_TYPE_DESC = Coursework:
contrast estimate SE df z.ratio p.value
Asian - Black 1.485 0.190 Inf 7.835 <.0001
Asian - Chinese 3.433 0.163 Inf 21.052 <.0001
Asian - Mixed -0.975 0.223 Inf -4.369 0.0001
Asian - White -2.288 0.127 Inf -18.035 <.0001
Black - Chinese 1.948 0.196 Inf 9.938 <.0001
Black - Mixed -2.459 0.247 Inf -9.971 <.0001
Black - White -3.773 0.165 Inf -22.835 <.0001
Chinese - Mixed -4.408 0.232 Inf -19.025 <.0001
Chinese - White -5.721 0.143 Inf -40.129 <.0001
Mixed - White -1.313 0.200 Inf -6.571 <.0001
ASSESSMENT_COMPONENT_TYPE_DESC = Exam:
contrast estimate SE df z.ratio p.value
Asian - Black 2.766 0.208 Inf 13.276 <.0001
Asian - Chinese 2.010 0.196 Inf 10.254 <.0001
Asian - Mixed -0.799 0.244 Inf -3.272 0.0094
Asian - White -2.626 0.139 Inf -18.861 <.0001
Black - Chinese -0.755 0.229 Inf -3.296 0.0087
Black - Mixed -3.564 0.270 Inf -13.203 <.0001
Black - White -5.391 0.181 Inf -29.788 <.0001
Chinese - Mixed -2.809 0.265 Inf -10.614 <.0001
Chinese - White -4.636 0.173 Inf -26.750 <.0001
Mixed - White -1.827 0.219 Inf -8.355 <.0001
Results are averaged over the levels of: PROGRAMME_PART
Degrees-of-freedom method: asymptotic
P value adjustment: tukey method for comparing a family of 5 estimates
#plot of mixed model 2
plot(mmodel2)
boxplot(mmodel2@resp$wtres~full_data2$ASSESSMENT_COMPONENT_TYPE_DESC)
boxplot(mmodel2@resp$wtres~full_data2$Ethnicity)
boxplot(mmodel2@resp$wtres~full_data2$PROGRAMME_PART)
Another way to check for homeoscedasticity - (sd) residuals
#ASSESSMENT TYPE
mmodel2res<-mmodel2@resp$wtres
#assessment type residuals in df
mmodel2df<-data.frame(full_data2$ASSESSMENT_COMPONENT_TYPE_DESC,mmodel2res)
names(mmodel2df)<-c("ASSESSMENT_COMPONENT_TYPE_DESC","residual")
mmodel2df %>% group_by(ASSESSMENT_COMPONENT_TYPE_DESC)%>% summarise(sd(residual))
The assumption is met because highest sd is less than twice the lowest.
#Ethnicity
mmodel2df2<-data.frame(full_data2$Ethnicity,mmodel2res)
names(mmodel2df2)<-c("Ethnicity","residual")
mmodel2df2 %>% group_by(Ethnicity)%>% summarise(sd(residual))
Again, the assumption is met.
#Programme part
mmodel2df3<-data.frame(full_data2$PROGRAMME_PART,mmodel2res)
names(mmodel2df3)<-c("Programme_Part","residual")
mmodel2df3 %>% group_by(Programme_Part)%>% summarise(sd(residual))
Assumption satisfied again..
Moving onto checking the residuals of the model are normally distributed.
#now check for residuals for mmodel2
qqnorm(mmodel2res)
#confint(mmodel2)
Data set for this question will be full_data data set as already includes Gender data. Created this in preliminary analysis.
Distribution of Marks of Male & Female students (2017-2021)
boxplot(MARK~Gender, data = full_data)
Below is bar chart showing of Gender of ethnicity on mean mark between 2017-2021
#Visualise effect of ethnicity on overall/mean mark between 2017-2021
full_data %>%
select(Gender, MARK) %>%
group_by(Gender) %>%
tally(mean(MARK)) %>%
rename(Average_mark = n) %>%
ggplot(aes(x=Gender, y = Average_mark, fill =Gender)) +
geom_bar(stat = 'identity', position = 'dodge') +
scale_fill_brewer(palette = "Paired") +
geom_text(aes(label = round(Average_mark, 2)), position=position_dodge(width=0.9), vjust=-0.25) +
ggtitle("Average Test Scores by Gender between 2017-2021") +
theme(plot.title = element_text(hjust = 0.5)) +
ylab('Average Mark')
#looking at avg grade for each ethnicity type between 2017-2021
full_data %>%
select(INSTANCE_ACADEMIC_YEAR, Gender,MARK) %>%
group_by(INSTANCE_ACADEMIC_YEAR, Gender) %>%
tally(mean(MARK)) %>%
rename(Average_mark = n) %>%
#spread(key = Ethnicity, value = Average_mark) %>%
ggplot(aes(x = INSTANCE_ACADEMIC_YEAR, y = Average_mark, col = Gender)) +
geom_line() +
geom_point() +
ggtitle("Average grade for Men vs Women (2017-2021)") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab('Academic Year') +
ylab('Average Mark') +
labs(col='Gender')
#Visualising percentage of white vs BAME students achieving top degree (2017-2021)
full_data %>%
group_by(ID2) %>%
mutate(Average_mark = mean(MARK)) %>%
filter(MARK >=60) %>%
select(ID2, Gender, Average_mark) %>%
distinct() %>%
group_by(Gender) %>%
summarise(counts = n()) %>%
mutate(percentage = scales::label_percent()(counts / sum(counts))) %>%
ggplot(aes(x=Gender, y = counts, fill =Gender)) +
geom_bar(stat = 'identity', position = 'dodge') +
scale_fill_brewer(palette = "Pastel1") +
geom_text(aes(label = percentage), position = position_stack(vjust = 1.1), size = 3) +
ggtitle("Number of Men vs Women achieving top degree (2017-2021)") +
theme(plot.title = element_text(hjust = 0.5)) +
ylab('Number of students')
This is very interesting as more men are being awarded top degrees but on average woman perform better. This may be because more men in this dataset then woman. So, below is a better comparison, looking at distribution of grades.
#Gender AWARDING GAP
full_data %>%
group_by(ID2) %>%
mutate(Average_mark = mean(MARK)) %>%
mutate(Grade = case_when(Average_mark >= 70 ~ "First Class",
Average_mark >= 60 ~ "Upper Second Class Honour",
Average_mark >= 50 ~ "Lower Second Class Honour",
Average_mark >= 40 ~ "Third class/pass",
Average_mark < 40 ~ "Fail")) %>%
group_by(Gender, Grade) %>%
summarise(counts = n()) %>%
mutate(freq = (counts / sum(counts))) %>%
mutate(percentage = scales::label_percent()(counts / sum(counts))) %>%
ggplot(aes(x=Gender, y = freq, fill =Grade)) +
geom_bar(stat = 'identity') +
coord_flip() +
scale_fill_brewer(palette = "Pastel1") +
geom_text(aes(label = percentage), position = position_stack(vjust = 0.5), size = 2.5) +
ggtitle("Distribution of degree outcomes by gender (2017-2021)") +
theme(plot.title = element_text(hjust = 0.5)) +
labs(y = "Frequency", fill = "Degree Classification")
`summarise()` has grouped output by 'Gender'. You can override using the `.groups` argument.
This graph is much better at showing what is actually going on in regards to what degree classification students are being awarded. This makes more sense now… This is showing 33.5% of woman students are getting a first in comparison to male students So females do perform better/ more females are awarded a top degree in comparison to male students, but when looking at overall degree outcomes, more males are getting top degree because there are more male students at university.
full_data %>%
select(Gender, ASSESSMENT_COMPONENT_TYPE_DESC,MARK) %>%
group_by(Gender, ASSESSMENT_COMPONENT_TYPE_DESC) %>%
tally(mean(MARK)) %>%
rename(Average_mark = n) %>%
ggplot(aes(x=ASSESSMENT_COMPONENT_TYPE_DESC, y = Average_mark, fill = Gender)) +
geom_bar(stat = 'identity', position = 'dodge') +
scale_fill_brewer(palette = "Set1") +
geom_text(aes(label = round(Average_mark, 2)), position=position_dodge(width=0.9), vjust=-0.25) +
ggtitle("Average Mark By Assessment Type And Gender") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab('Assessment Type') +
ylab('Average Mark')
This shows women perform better than men in both assessment types, but the difference by assessment type seems similar… So one gender does not perform significantly different than other in one assessment type.
full_data %>%
group_by(ID2) %>%
mutate(Average_mark = mean(MARK)) %>%
filter(Average_mark >=60) %>%
select(ID2, Gender, ASSESSMENT_COMPONENT_TYPE_DESC, Average_mark) %>%
distinct() %>%
select(Gender, ASSESSMENT_COMPONENT_TYPE_DESC) %>%
group_by(Gender,ASSESSMENT_COMPONENT_TYPE_DESC) %>%
summarise(counts = n()) %>%
mutate(percentage = scales::label_percent()(counts / sum(counts))) %>%
#select(Gender, ASSESSMENT_COMPONENT_TYPE_DESC, percentage) %>%
ggplot(aes(x=ASSESSMENT_COMPONENT_TYPE_DESC, y = percentage, fill = Gender)) +
geom_bar(stat = 'identity', position = 'dodge') +
scale_fill_brewer(palette = "Set1") +
geom_text(aes(label = percentage), position=position_dodge(width=0.9), vjust=-0.25) +
ggtitle("Effect of assessment type on gender awarding gap") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab('Assessment Type') +
ylab('Percentage of students')
Adding missing grouping variables: `ID2`
`summarise()` has grouped output by 'Gender'. You can override using the `.groups` argument.
5% diff in each assessment type… Doesnt look like exams are leading to a wider gender awarding gap. Difference in performance for exams is same as difference for coursework….
interaction.plot(x.factor = full_data$Gender, #x-axis variable
trace.factor = full_data$ASSESSMENT_COMPONENT_TYPE_DESC, #variable for lines
response = full_data$MARK, #y-axis variable
fun = median, #metric to plot
ylab = "Mark",
xlab = "Ethnicity",
col = c("green", "blue"),
lty = 1, #line type
lwd = 2, #line width
trace.label = "Assessment type")
# **Research Question 3**
# **3. Is there a link between assessment type, gender and student achievement?**
mmodel3<- lmer(MARK ~ ASSESSMENT_COMPONENT_TYPE_DESC * Gender + PROGRAMME_PART + (1|ID2) + (1|MODULE_CODE), data = full_data)
summary(mmodel3)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: MARK ~ ASSESSMENT_COMPONENT_TYPE_DESC * Gender + PROGRAMME_PART +
(1 | ID2) + (1 | MODULE_CODE)
Data: full_data
REML criterion at convergence: 6689154
Scaled residuals:
Min 1Q Median 3Q Max
-7.4574 -0.5327 0.0223 0.5787 4.7917
Random effects:
Groups Name Variance Std.Dev.
ID2 (Intercept) 40.57 6.369
MODULE_CODE (Intercept) 25.44 5.044
Residual 153.51 12.390
Number of obs: 840368, groups: ID2, 31540; MODULE_CODE, 6945
Fixed effects:
Estimate Std. Error df t value
(Intercept) 6.729e+01 1.431e-01 1.236e+04 470.102
ASSESSMENT_COMPONENT_TYPE_DESCExam -8.194e+00 4.717e-02 7.750e+05 -173.732
GenderWoman 8.882e-01 8.896e-02 3.505e+04 9.985
PROGRAMME_PARTB 2.593e-02 1.613e-01 1.548e+04 0.161
PROGRAMME_PARTC 7.156e-01 1.899e-01 1.292e+04 3.767
PROGRAMME_PARTD -1.454e+00 3.780e-01 3.142e+04 -3.846
PROGRAMME_PARTF 5.782e+00 5.344e-01 6.410e+03 10.819
PROGRAMME_PARTT -1.546e+00 1.887e-01 1.131e+04 -8.194
ASSESSMENT_COMPONENT_TYPE_DESCExam:GenderWoman 8.917e-01 6.904e-02 8.360e+05 12.916
Pr(>|t|)
(Intercept) < 2e-16 ***
ASSESSMENT_COMPONENT_TYPE_DESCExam < 2e-16 ***
GenderWoman < 2e-16 ***
PROGRAMME_PARTB 0.872332
PROGRAMME_PARTC 0.000166 ***
PROGRAMME_PARTD 0.000120 ***
PROGRAMME_PARTF < 2e-16 ***
PROGRAMME_PARTT 2.8e-16 ***
ASSESSMENT_COMPONENT_TYPE_DESCExam:GenderWoman < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) ASSESSMENT_COMPONENT_TYPE_DESCEx GndrWm
ASSESSMENT_COMPONENT_TYPE_DESCEx -0.117
GenderWoman -0.245 0.104
PROGRAMME_PARTB -0.580 -0.005 -0.003
PROGRAMME_PARTC -0.557 0.000 0.008
PROGRAMME_PARTD -0.297 0.006 0.014
PROGRAMME_PARTF -0.230 0.005 0.014
PROGRAMME_PARTT -0.696 0.036 -0.023
ASSESSMENT_COMPONENT_TYPE_DESCE: 0.048 -0.467 -0.209
PROGRAMME_PARTB PROGRAMME_PARTC PROGRAMME_PARTD
ASSESSMENT_COMPONENT_TYPE_DESCEx
GenderWoman
PROGRAMME_PARTB
PROGRAMME_PARTC 0.510
PROGRAMME_PARTD 0.246 0.328
PROGRAMME_PARTF 0.153 0.144 0.075
PROGRAMME_PARTT 0.444 0.430 0.255
ASSESSMENT_COMPONENT_TYPE_DESCE: 0.003 0.000 -0.001
PROGRAMME_PARTF PROGRAMME_PARTT
ASSESSMENT_COMPONENT_TYPE_DESCEx
GenderWoman
PROGRAMME_PARTB
PROGRAMME_PARTC
PROGRAMME_PARTD
PROGRAMME_PARTF
PROGRAMME_PARTT 0.170
ASSESSMENT_COMPONENT_TYPE_DESCE: -0.002 0.003
This is interesting…shows gap in AT smaller for woman compared to men. Also average cw grade very similar for men and woman And woman perform better in exams compared to men And woman perform better in exams compared to men…
plot_model(mmodel3, type = "pred", terms = c("Gender", "ASSESSMENT_COMPONENT_TYPE_DESC"))
anova(mmodel3)
Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
ASSESSMENT_COMPONENT_TYPE_DESC 4858953 4858953 1 716363 31652.242 < 2.2e-16 ***
Gender 34941 34941 1 34071 227.611 < 2.2e-16 ***
PROGRAMME_PART 44305 8861 5 12859 57.723 < 2.2e-16 ***
ASSESSMENT_COMPONENT_TYPE_DESC:Gender 25610 25610 1 835965 166.826 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Again, similar to RQ1 & RQ2, all variables hold a significant effect on student performance
#effect size...look at video for how to explain this in report.
effectsize::eta_squared(mmodel3, partial = TRUE)
# Effect Size for ANOVA (Type III)
Parameter | Eta2 (partial) | 95% CI
---------------------------------------------------------------------
ASSESSMENT_COMPONENT_TYPE_DESC | 0.04 | [0.05, 1.00]
Gender | 6.64e-03 | [0.01, 1.00]
PROGRAMME_PART | 0.02 | [0.02, 1.00]
ASSESSMENT_COMPONENT_TYPE_DESC:Gender | 2.00e-04 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
#emmeans is meant to be looked at for main effect.
#in this case they were all significant so will do for each??
emmeans(mmodel3, list(pairwise~ASSESSMENT_COMPONENT_TYPE_DESC), adjust="tukey")
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'pbkrtest.limit = 840368' (or larger)
[or, globally, 'set emm_options(pbkrtest.limit = 840368)' or larger];
but be warned that this may result in large computation time and memory use.
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'lmerTest.limit = 840368' (or larger)
[or, globally, 'set emm_options(lmerTest.limit = 840368)' or larger];
but be warned that this may result in large computation time and memory use.
NOTE: Results may be misleading due to involvement in interactions
$`emmeans of ASSESSMENT_COMPONENT_TYPE_DESC`
ASSESSMENT_COMPONENT_TYPE_DESC emmean SE df asymp.LCL asymp.UCL
Coursework 68.32 0.1268 Inf 68.07 68.57
Exam 60.57 0.1297 Inf 60.32 60.83
Results are averaged over the levels of: Gender, PROGRAMME_PART
Degrees-of-freedom method: asymptotic
Confidence level used: 0.95
$`pairwise differences of ASSESSMENT_COMPONENT_TYPE_DESC`
1 estimate SE df z.ratio p.value
Coursework - Exam 7.75 0.0436 Inf 177.911 <.0001
Results are averaged over the levels of: Gender, PROGRAMME_PART
Degrees-of-freedom method: asymptotic
emmeans(mmodel3, list(pairwise~Gender), adjust="tukey")
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'pbkrtest.limit = 840368' (or larger)
[or, globally, 'set emm_options(pbkrtest.limit = 840368)' or larger];
but be warned that this may result in large computation time and memory use.
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'lmerTest.limit = 840368' (or larger)
[or, globally, 'set emm_options(lmerTest.limit = 840368)' or larger];
but be warned that this may result in large computation time and memory use.
NOTE: Results may be misleading due to involvement in interactions
$`emmeans of Gender`
Gender emmean SE df asymp.LCL asymp.UCL
Man 63.78 0.1301 Inf 63.53 64.04
Woman 65.12 0.1376 Inf 64.85 65.38
Results are averaged over the levels of: ASSESSMENT_COMPONENT_TYPE_DESC, PROGRAMME_PART
Degrees-of-freedom method: asymptotic
Confidence level used: 0.95
$`pairwise differences of Gender`
1 estimate SE df z.ratio p.value
Man - Woman -1.33 0.0884 Inf -15.087 <.0001
Results are averaged over the levels of: ASSESSMENT_COMPONENT_TYPE_DESC, PROGRAMME_PART
Degrees-of-freedom method: asymptotic
emmeans(mmodel3, list(pairwise~PROGRAMME_PART), adjust="tukey")
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'pbkrtest.limit = 840368' (or larger)
[or, globally, 'set emm_options(pbkrtest.limit = 840368)' or larger];
but be warned that this may result in large computation time and memory use.
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'lmerTest.limit = 840368' (or larger)
[or, globally, 'set emm_options(lmerTest.limit = 840368)' or larger];
but be warned that this may result in large computation time and memory use.
$`emmeans of PROGRAMME_PART`
PROGRAMME_PART emmean SE df asymp.LCL asymp.UCL
A 63.86 0.1384 Inf 63.59 64.13
B 63.89 0.1354 Inf 63.62 64.15
C 64.58 0.1584 Inf 64.27 64.89
D 62.41 0.3612 Inf 61.70 63.12
F 69.64 0.5200 Inf 68.62 70.66
T 62.31 0.1309 Inf 62.06 62.57
Results are averaged over the levels of: ASSESSMENT_COMPONENT_TYPE_DESC, Gender
Degrees-of-freedom method: asymptotic
Confidence level used: 0.95
$`pairwise differences of PROGRAMME_PART`
1 estimate SE df z.ratio p.value
A - B -0.0259 0.161 Inf -0.161 1.0000
A - C -0.7155 0.190 Inf -3.767 0.0023
A - D 1.4538 0.378 Inf 3.846 0.0017
A - F -5.7823 0.534 Inf -10.819 <.0001
A - T 1.5463 0.189 Inf 8.194 <.0001
B - C -0.6896 0.176 Inf -3.927 0.0012
B - D 1.4797 0.373 Inf 3.969 0.0010
B - F -5.7563 0.534 Inf -10.777 <.0001
B - T 1.5723 0.186 Inf 8.450 <.0001
C - D 2.1693 0.363 Inf 5.974 <.0001
C - F -5.0667 0.541 Inf -9.370 <.0001
C - T 2.2619 0.202 Inf 11.187 <.0001
D - F -7.2361 0.631 Inf -11.469 <.0001
D - T 0.0925 0.377 Inf 0.245 0.9999
F - T 7.3286 0.536 Inf 13.681 <.0001
Results are averaged over the levels of: ASSESSMENT_COMPONENT_TYPE_DESC, Gender
Degrees-of-freedom method: asymptotic
P value adjustment: tukey method for comparing a family of 6 estimates
#There was an interaction so need emmeans for this
emmeans(mmodel3, pairwise ~ Gender | ASSESSMENT_COMPONENT_TYPE_DESC)
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'pbkrtest.limit = 840368' (or larger)
[or, globally, 'set emm_options(pbkrtest.limit = 840368)' or larger];
but be warned that this may result in large computation time and memory use.
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'lmerTest.limit = 840368' (or larger)
[or, globally, 'set emm_options(lmerTest.limit = 840368)' or larger];
but be warned that this may result in large computation time and memory use.
$emmeans
ASSESSMENT_COMPONENT_TYPE_DESC = Coursework:
Gender emmean SE df asymp.LCL asymp.UCL
Man 67.88 0.1307 Inf 67.62 68.13
Woman 68.77 0.1380 Inf 68.50 69.04
ASSESSMENT_COMPONENT_TYPE_DESC = Exam:
Gender emmean SE df asymp.LCL asymp.UCL
Man 59.68 0.1337 Inf 59.42 59.95
Woman 61.46 0.1442 Inf 61.18 61.75
Results are averaged over the levels of: PROGRAMME_PART
Degrees-of-freedom method: asymptotic
Confidence level used: 0.95
$contrasts
ASSESSMENT_COMPONENT_TYPE_DESC = Coursework:
contrast estimate SE df z.ratio p.value
Man - Woman -0.888 0.089 Inf -9.985 <.0001
ASSESSMENT_COMPONENT_TYPE_DESC = Exam:
contrast estimate SE df z.ratio p.value
Man - Woman -1.780 0.101 Inf -17.704 <.0001
Results are averaged over the levels of: PROGRAMME_PART
Degrees-of-freedom method: asymptotic
Do not need to check linearity as each variable is categorical.
#plot of mixed model 3 residuals
plot(mmodel3)
boxplot(mmodel3@resp$wtres~full_data$ASSESSMENT_COMPONENT_TYPE_DESC)
boxplot(mmodel3@resp$wtres~full_data$Gender)
boxplot(mmodel3@resp$wtres~full_data$PROGRAMME_PART)
Another way to check for homeoscedasticity - (sd) residuals
#ASSESSMENT TYPE
mmodel3res<-mmodel3@resp$wtres
#assessment type residuals in df
mmodel3df<-data.frame(full_data$ASSESSMENT_COMPONENT_TYPE_DESC,mmodel3res)
names(mmodel3df)<-c("ASSESSMENT_COMPONENT_TYPE_DESC","residual")
mmodel3df %>% group_by(ASSESSMENT_COMPONENT_TYPE_DESC)%>% summarise(sd(residual))
#Gender
mmodel3df2<-data.frame(full_data$Gender,mmodel3res)
names(mmodel3df2)<-c("Gender","residual")
mmodel3df2 %>% group_by(Gender)%>% summarise(sd(residual))
#Programme part
mmodel3df3<-data.frame(full_data$PROGRAMME_PART,mmodel3res)
names(mmodel3df3)<-c("Programme_Part","residual")
mmodel3df3 %>% group_by(Programme_Part)%>% summarise(sd(residual))
Assumptions met for all variables
qqnorm(mmodel3res)
All assumptions for mixed model 3 met.
#confint(mmodel3)